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ABSTRACT 


A  hybrid-stress  formulation  of  isoparametric  elements  for  the  analysis 
of  thin  multilayer  laminated  composite  plates  is  presented,  and  is  applied  to 
edge  contact  analyses.  The  element  displacement  behavior  is  characterized 
by  laminate  reference  surface  inplane  and  transverse  displacements  and 
laminate  non-normal  cross-section  rotations;  as  a  result,  the  number  of 
degrees  of  freedom  is  independent  of  the  number  of  layers.  All  components 
of  stress  are  included  and  are  related  to  a  set  of  laminate  stress  param¬ 
eters,  the  number  of  which  is  independent  of  the  number  of  layers.  Atten¬ 
tion  is  restricted  here  to  thin  laminates:  for  thin  laminates  it  is  shown 
that  the  contributions  of  transverse  shear  stress  and  transverse  normal 
stress  to  the  internal  complementary  strain  energy  £an  be  neglected.  As  a 
result,  a  modified  stiffness-formation-algorithm  can  be  used  which  provides 
a  significant  improvement  in  computation  efficiency.  The  formulation  is  used 
to  develop  an  8-node  isoparametric  thin  multilayer  plate  element.  The 
resulting  element  is  naturally  invariant,  of  correct  rank,  and  non-locking 
in  the  thin  plate  limit.  Element  performance  is  documented  here  jjor  several 
illustrative  examples.  / ,  ,  ,  •  /  /—  -  -r  ,  ,  •  /  r  , 

The  newly  developed  Element  is  applied  to  the  analysis  of  laminate  edge 
contact.  A  hybrid-stress  formulation  for  the  incremental/iterative  analysis 
of  contact  problems  as  well  as  the  procedure  to  locate  the  surface  of  contact 
is  presented.  The  assumed  contact  surface  is  divided  into  contact  elements 
having  unknown  nodal  contact  tractions.  The  finite  element  equations  are 
solved  for  the  element  nodal  degrees  of  freedom  as  well  as  unknown  nodal 
contact  tractions.  The  assumed  contact  surface  is  checked  to  determine  if  it 
satisfies  the  contact  conditions.  If  not  a  new  contact  surface  is  assumed 
and  the  (iterative)  process  is  repeated.  Attention  is  restricted  here  to 
symmetric  laminates  subject  to  contact  in  directions  normal  to  the  laminate 
edge,  so  that  only  inplane  displacement  occurs.  Several  examples  involving 
elastic,  non- frictional  contracting  bodies  are  presented  to  verify  the 
formulation/algorithm.  An  example  involving  a  laminated  plate  is  also 
presented. 
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Chapter  I 
INTRODUCTION 

Multilayer  £iber  reinforced  conposites  are  of  interest  in  a  number 
of  structural  applications  where  high  strength  to  weight  ratios  are  re¬ 
quired.  An  area  of  recent  and  continuing  interest  is  the  failure  analy¬ 
sis  of  composite  laminates.  Composites  are  strongly  anisotropic  by  na- 

m 

ture  and  contain  major  microscopic  inhomogeneities  corresponding  to  the 
matrix/fiber  nature  of  the  materials.  Therefore  accurate  and  detailed 
stress  distributions  must  be  available  in  order  to  predict  failure  of 
laminated  composite  structures. 

Most  engineering  devices  are  made  up  of  assembled  parts  and  many  of 
these  are  mechanically  joined.  Even  under  small  deformations  nonlinear¬ 
ities  often  arise  in  these  joints.  The  nonlinearities  arise  for  in¬ 
stance  in  contact  problems  such  as  bolted  connections,  shrink  fit,  or 
roller  bearings.  Sometimes  gaps  separate  various  structural  components 
presenting  nonlinear  problems  which  involve  rigid  body  movement  across  a 
gap  with  subsequent  contact  of  disconnected  regions.  In  many  situations 
structural  components  in  contact  slide  relative  to  one  another. 

The  complicating  factor  in  these  problems  is  the  unknown  surface  of 
contact.  In  aany  cases,  it  is  important  to  know  the  contact  conditions 
in  order  to  predict  accurately  the  stresses  and  strength  of  the  joint. 

Analytic  solutions  of  contact  problems  are  limited  to  idealized 
simple  configurations,  loadings,  and  boundary  conditions.  These  solu- 
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tions  may  be  applied  to  design  problems  with  varying  degrees  of  success 
depending  upon  how  closely  the  geometry  and  loading' model  the  actual 
problem.  Thus,  numerical  techniques  are  needed  for  solutions  of  realis¬ 
tic  contact  problems. 

The  finite  element  method  will  be  used  to  develop  an  isoparametric 
multilayer  plate  element  based  on  a  hybrid-stress  formulation.  The  dis¬ 
placement  behavior  will  be  characterized  by  laminate  reference  plane  in¬ 
plane  and  transverse  displacements  and  laminate  non-normal  cross-section 
rotations;  as  a  result,  the  number  of  degrees  of  freedom  will  be  inde¬ 
pendent  of  the  number  of  layers.  All  components  of  stress  will  be  in¬ 
cluded. 

This  multilayer  plate  element  will  be  used  to  solve  edge  contact 
problems.  A  hybrid-stress  formulation  for  the  analysis  of  contact  prob¬ 
lems  as  well  as  the  procedure  to  locate  the  surface  of  contact  will  be 
presented.  The  contact  problem  may  be  either  frictional  or  frictionless 
and  may  involve  extensive  sliding  between  deformable  bodies. 

Chapter  II  is  devoted  to  the  development  of  the  multilayer  plate 
element.  A  review  of  recent  work  is  included  in  section  2.1.  A  stiff¬ 
ness  matrix  is  formulated  for  a  moderately  thick  multilayer  plate  ele¬ 
ment  in  section  2.2.  This  stiffness  is  reduced  to  obtain  a  stiffness 
matrix  for  thin  multilayer  plate  elements  in  section  2.3.  Section  2.4 
describes  the  formation  and  evaluation  of  a  8-node  thin  multilayer  plate 
element. 

Chapter  III  describes  the  hybrid-stress  formulation  for  the  analy¬ 
sis  of  contact  problems  as  well  as  the  solution  technique  and  applica¬ 
tions.  A  review  of  recent  work  is  included  in  section  3.1.  The  hybrid 
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streat  functional  for  contact  it  presented  in  section  3.2.  The  matrix 
equations  are  given  in  section  3.3.  Section  3.4  describes  the  contact 
eleaient  and  section  3.5  develops  the  solution  technique.  Application 
and  evaluation  are  included  in  section  3.6. 

Chapter  IV  i*  a  statement  of  the  summary  and  conclusions  of  this 
work.  This  chapter  concludes  with  some  suggestions  for  further  work. 
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Chapter  II 

HYBRID  STRESS  REDUCED-MINDLIN  ELEMENTS  FOR  THIN  MULTILAYER 

PLATES 


2.1  INTRODUCTION 

Multilayer  fiber  reinforced  compoaites  are  of  intereat  in  a  number 
of  atructural  applicationa  where  high  atrength  to  weight  ratios  are  re¬ 
quired.  A  number  of  multilayer  plate  elements  have  appeared  in  the  lit¬ 
erature,  ranging  in  applicability  from  thin  (lamination  theory)  to  mod¬ 
erately  thick  (including  transverse  shear  effects)  laminated  plates 
(e.g.  Cl— 1 1] ) .  These  elements  are  based  prirarily  on  assumed-displace¬ 
ment  or  hybrid-  stress  formulations,  although  elements  based  on  other 
variational  principles  have  been  proposed  (e.g.  [12-13]). 

There  are  advantages  to  the  use  of  the  hybrid-stress  model  for  mul¬ 
tilayer  plate  elements.  Hybrid-stress  elements  can  be  derived  which  in¬ 
clude  individual  (independent)  layer  cross-section  rotations  (also  pos¬ 
sible  in  assumed-displacement  elements)  and  individual (independent) 
layer  stress  fields.  The  layer  stress  fields  can  be  selected  such  that 
interlayer  surface  traction  continuity  and  laminate  upper/lower  surface- 
traction  free  conditions  are  exactly  satisfied;  examples  of  hybrid- 
stress  elements  of  this  type  are  found  in  References  [6,10,11].  Ele¬ 
ments  in  this  category  are  applicable  to  moderately  thick  and  thick 
laminates  (i.e.  typical  thickness  to  span  ratios  as  high  as  h/L-0.25). 
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Second,  by  proper  choice  of  the  spanwise  distributions  of  the  in¬ 
plane  stresses,  '  it  is  possible  to  derive  single-layer  hybrid-stress 
plate  elements  which  are  based  on  Mindlin-type  through  thickness  dis¬ 
placements  (i.e.  independent  non-normal  cross-section  rotations)  and  in¬ 
clude  all  components  of  stress,  and  which  are  of  correct  rank  and  non¬ 
locking  in  the  thin-plate  limit  [14-16] .  ("Locking"  refers  here  to 
excessively  stiff  solutions  obtained  with  some  moderately  thick  plate 
elements  when  applied  to  thin  plate  problems.)  The  characteristic  of 
natural  invariance  can  also  be  achieved  for  the  8-node  plate  element 
[17].  These  attributes  carry  over  to  hybrid-stress  multilayer  plate  el¬ 
ements  [10,11];  the  8-node  element  of  Reference  [ll] ,  for  example,  is 
accurate  for  moderately  thick  to  thick  laminates,  is  of  correct  rank  and 
naturally  invariant,  and  is  non-locking  when  used  for  thin  laminated 
plates. 

In  many  applications  of  laminated  composites,  the  plate  can  be  as¬ 
sumed  to  be  thin  and  thus  governed  by  lamination  theory.  Application  of 
elements  of  the  type  in  References  [10,11],  for  example,  to  thin  plates 
may  be  considered  inefficient.  In  those  elements  the  number  of  degrees 
of  freedom  and  the  number  of  stress  parameters  will  grow  in  proportion 
to  the  number  of  layers,  resulting  in  increased  computation  times  and 
storage  requirements  for  element  stiffness  formation.  This  generality 
at  the  expense  of  significantly  increased  computation  cost  is  therefore 
not  warranted.  One  approach  used  in  References  [10,11]  was  to  kinemati¬ 
cally  constrain  the  through-thickness  displacement  field  so  that  element 
nodal  degrees  of  freedom  corresponded  to  reference  plane  displacements 
and  laminate  cross-section  rotations  -  independent  of  the  number  of 


6 


layers.  The  eleaents  produced  accurate  predictions  for  typical  thick¬ 
ness  ratios  h/L  -0.1  and  lower.  Those  "restricted"  eleaents  [10,11] 
continued  to  be  based  on  independent  layer  stress  fields  and  thus  eoapu- 
tation  tiae  and  storage  requireaents  for  the  foraation  of  an  eleaent 
stiffness  continued  to  grow  with  the  nuaber  of  layers.  Furtheraore,  the 
use  of  independent  layer  stress  fields  with  a  laainate  displaceaent 
field  seeas  inconsistent. 

The  foraulation  presented  here  rectifies  this  inconsistency  by  pro¬ 
ducing  hybrid-stress  aultilayer  plate  eleaents  for  which  the  nuaber  of 
nodal  degrees  of  freedoa  and  the  nuaber  of  stress  paraaeters  is  indepen¬ 
dent  of  the  nuaber  of  layers.  In  these  eleaents,  degrees  of  freedoa  at 
a  node  correspond  to  laainate  reference  surface  inplane  and  transverse 
displaceaents  and  laainate  cross-section  rotations  (i.e.  a  total  of  5 
dof  per  node  independent  of  the  nuaber  of  layers) .  In  order  to  produce 
stress  fields  in  each  layer  related  to  a  fixed  set  of  laainate  stress 
paraaeters,  laainate  inplane  strains  are  first  interpolated  in  terns  of 
stress  paraaeters.  Layer  inplane  stress  fields  are  related  to  these  pa¬ 
raaeters  via  constitutive  equations.  Layer  transverse  shear  and  tran¬ 
sverse  noraal  stresses  are  obtained  froa  the  inplane  stesses  via  inte¬ 
gration  of  the  hoaogeneous  equilibriua  equations;  the  constants  of 
integration  are  selected  to  satisfy  lower  surface  traction-free  condi¬ 
tions  and  interlayer  surface  traction  continuity. 

This  aproach  was,  in  effect,  used  for  the  developaent  of  a  4-node 
hybrid-stress  aultilayer  eleaent  in  Reference  [9] .  Results  in  that 
study  suggest  that  aoderately-thick  (i.e.  h/L20.1)  laainates  can  be  ac¬ 


curately  analyzed 


The  mu'*'ilayer  element  behavior,  particularly  in  the  thin  plate 
limit,  will  depend  on  the  form  of  the  spanwise  distribution  of  the 
stresses,  just  as  in  the  single  layer  hybrid-stress  Mindlin  type  ele¬ 
ments.  Thus,  early  phases  of  the  present  work  aimed  at  the  development 
of  an  8-node  multilayer  element  (with  laminate  degrees  of  freedom  and 
stress  parameters)  having  a  spanwise  stress  distribution  analogous  to 
the  single  layer  element  QH3  of  Reference  [17].  The  resulting  multilay¬ 
er  element  was  of  correct  rank,  non-locking,  and  naturally  invariant. 
However  computation  times  became  excessive  (as  will  be  described  later) 
when  using  the  approach  of  Reference  [9]  for  large  numbers  of  layers; 
this  motivated  the  search  for  an  alternate  approach  for  defining  a  hy¬ 
brid-stress  multilayer  plate  element. 

In  a  recent  study  [18],  Spilker  has  shown  that  the  single-layer  hy¬ 
brid-stress  elements  of  References  [14-17]  can  be  reduced  to  thin  plate 
elements  by  neglecting  the  contributions  of  the  transverse  shear  and 
transverse  normal  stresses  to  the  internal  complementary  energy.  Sig¬ 
nificant  computational  savings  were  found  while  maintaining  the  ease  of 
formulating  plate  elements  based  on  independent  transverse  displacement 
and  cross-section  rotations. 

The  present  study  extends  the  concepts  of  Reference  [18]  to  multi¬ 
layer  plates.  Described  herein  is  a  hybrid-stress  formulation  for  thin 
multilayer  laminated  plates  with  nodal  degrees  of  freedom  corresponding 
to  laminate  reference  surface  inplane  and  transverse  displacements  and 
laminate  cross-section  rotations.  All  components  of  stress  are  included 
and  are  related  to  a  set  of  laminate  stress  parameters,  with  contribu¬ 
tions  of  transverse  shear  and  normal  stresses  to  the  internal  complemen- 
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tary  energy  neglected.  Computation  time  for  the  resulting  element 
stiffness  aatrix  is  shown  to  be  essentially  independent  of  the  number  of 
layers. 

To  illustrate  the  foraulation,  an  8-node  multilayer  plate  element 
is  presented  which  is  analogous  to  the  single  layer  element  QH3T  [18] . 
The  8-node  element  is  naturally  invariant,  of  correct  rank,  and  non¬ 
locking.  Results  for  several  example  problems  show  that  the  element 
converges  to  the  thin  plate  (lamination  theory)  solution. 


2.2  STIFFNESS  MATRIX  FOR  A  MODERATELY  THICK  MULTILAYER  PLATE  ELEMENT 

The  multilayer  plate  is  assumed  to  lie  in  the  x-y  plane,  with  z  ■  0 
corresponding  to  the  laminate  reference  surface,  located  arbitrarily 
(e.g.  geometric  midsurface).  The  laminate  consists  of  I  perfectly-bond¬ 
ed  layers  numbered  bottom  to  top,  with  z  *  h,»  h2  ...  hj+j  locat¬ 
ing  the  lower,  interlayer,  and  upper  surfaces,  respectively  (see  Figure 
la) .  Displacement  behavior  is  characterized  by  laminate  reference  sur¬ 
face  displacements  uq,  vq,  and  wq,  in  the  x,y,  and  z  directions, 
respectively,  and  laminate,  cross-section  rotations,  0X  and  6y 
(which  are  assumed  to  be  independent  of  w) ,  about  the  x  and  y  axes,  re¬ 
spectively  (see  Figure  lb).  Note  that  inplane  displacements  are  includ¬ 
ed  due  to  inherent  material  induced  bending/stretching  coupling  in  mul¬ 
tilayer  laminated  plates. 

The  hybrid-stress  functional  for  multilayer  plates,  assuming  per¬ 
fectly  bonded  layers  and  traction  continuity  on  interlayer  surfaces,  is 
given  by  [10] . 


(a)  layer  nuabering,  interface  coordinates 


Figure  1  Nodal  degrees  of  freedoafor  the 
through'thicknesa  displaceaent  field. 
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7T*.  zz[k\  s‘rf£'  JV  -  j  "2  dv 

(2.1) 

+■  ^  vjpclA 

A  »i 

where  aT- [crx  ay  az  ax2  gtg  ayz]  art  components  of  stress, 

£**[«x  ty  <z  «xy  <xz  tyj]  are  components  of  "strain"  as 
computed  from  displacements  via  the  linear  strain-  displacement  rela¬ 
tions.  Subscripts/superseripts  i  and  n  refer  to  layer  i  and  element  n, 
so  that  Vn£  is  the  volume  of  layer  i  for  element  n.  The  jS1  matrix 
relates  stresses  and  strains  in  the  xyz  coordinate  system  and  is  given 
in  detail  in  Reference  [10] .  The  last  term  in  equation  (2.1)  is  a  load 
term  corresponding  to  transverse  distributed  load,  p(x,y),  acting  over 
the  spanwise  area,  Ag,  of  element  n. 

The  purpose  of  the  present  section  is  to  define  the  form  of  the 
displacement  and  stress  interpolations  for  a  thin-to-moderately-thick 
multilayer  plate  element  for  which  the  number  of  nodal  degrees  of  free¬ 
dom  and  the  number  of  stress  parameters  is  independent  of  the  number  of 
layers.  Intermediate  element  matrices  and  the  element  stiffness  matrix 
will  also  be  defined.  However,  the  present  formulation  serves  as  the 
basis  for  the  development  of  efficient  thin  multilayer  plate  elements 
via  a  reduction  scheme  analogous  to  Reference  [18] .  The  reduction 
scheme  and  corresponding  elements  are  discussed  in  subsequent  subsec¬ 
tions. 


The  displacements  u,  v,  and  w  at  any  x,y,z  location  of  the  laminate 


are  given  by 


Ulx'h)l 


A) 


*4  <■*>  j)  -  y  (*>  ‘i ) 

£*,  i)  - 


(2.2) 


Using  the  linear  strain-displacement  relations,  the  inplane  strains. 


«Bl,  for  layer  i  are  given  by 

me  F 


A  i 

fa 

^  i 

f  iiUl 

/ax 

1 

’  Y  2  J 

) 

4  *  | 

. 

A*  1 

1  H  X 

1  ‘ 

X!f! 

!  a*u/  0  <?-o.y 

.  +  «* 

u  | 

h/«-^  i 

(2.3) 


The  transverse  shear  strains,  $j.*,  *r*  gi  v#n  by 


a  1 


4 


A  | 

Yx*, 

A  , 

[ 

*4 

^0, 


A*  + 


(2.  A) 


and  the  transverse  normal  strain,  $9X,  is  zero.  It  is  understood 
that  strains  in  layer  i  for  equations  (2.3)  and  (2. A)  are  obtained  by 
evaluating  z  between  h^  and  h^. 
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Denote  the  nodal  degrees  of  freedom  corresponding  to  uo  by 
thoaa  corresponding  to  vo  by^vo,  thosa  corresponding  to  wo  byj»o,  etc. 
The  spanwise  distributions  of  uo,  vo,  wo,  0X,  and  8y  can  be  interpo¬ 
lated  in  terms  of  nodal  degrees  of  freedom  by  using  any  of  the  faailies 
of  2-D  C*  continuity  shape  functions.  As  the  present  aultilayer  ele¬ 
ments  will  be  of  isoparametric  planform,  the  x,y  coordinates  are  related 
to  coordinates  |,»?  by  the  relation 

2  U<  CK,' J) 

i  '  ‘  (2.5) 

where  the  sum  extends  over  element  nodes,  (xj,  yj)  are  the  coordi¬ 
nates  of  element  node  j,  and  the  correspond  to  the  appropriate 

Serendipity  shape  functions. 

He  assume  here  that  the  shape  functions  used  for  the  displacements 
and  rotations  are  the  same  and  are  indentical  to  those  used  in  the  geo¬ 
metric  mapping  of  equation  (2.5).  Then 


represent  the  displacement/rotation  interpolations,  with  uoj,  v o  , 


etc.  corresponding  to  degrees  of  freedom  at  node  j.  For  convenience,  we 
define  the  interpolation  vector  N  and  element  nodal  degrees  of  freedom, 
q.  as 


(2.7) 


r=  L  uT  u.r  *jaT  exT 

-t)  ~  ^ 


Substituting  equation  (2.6)  into  equations  (2.3)  and  (2.4)  produces 


or 


'L> 

0 

6? 

o 

( - 

A 

% 

o 

n 

7 

0 

0 

o 

0 

D 

0 

fj 

?■» 

o 

0 

ox 

0 

4lfce 

-*% 

o 

t  -*£/ 

>  /a3 

Kl 

0 

JfcJ/ 

% 

-  Ki 

•V 

o 

— 

Is  (2.8) 


S'  -  ~  * 


(2.9) 


Note  that  N  “  N($,ij)  and  thus  derivatives  in  equation  (2.8)  will  incor- 
porate  the  mapping  of  equation  (2.5)  as  in  all  isoparametric  elements. 


The  stresses  gj-  in  each  layer  are  to  be  expressed  in  teras  of  a 
set  of  laainate  stress  paraaeters,  3,  the  total  nuaber  ef  which  is  inde¬ 
pendent  of  the  nuaber  of  layers  (i.e.  consistent  with  a  aultilayer  ele- 
aent  having  laainate  degrees  of  freedoa  independent  of  the  nuaber  of 
layers).  This  can  be  aecoaplished  by  using  a  scheme  siailar  to  that 
used  by  Spilker  et  al  [9]  for  a  4-node  hybrid-stress  aultilayer  ele- 
aent. 

The  laainate  inplane  strains  are  assuaed  to  vary  linearly  through 
the  laainate  thickness  and  are  expressed  in  teras  of  laainate  "stress 
paraaeters",  3,  in  the  fora 

sZk-  :i?b  J  ,  &  l  o.io) 

fcl 

In  equation  (2.10),  the  subscripts  s  and  b  refer  to  stretching  and  bend¬ 
ing  contributions,  respectively.  The  interpolation  matrices  7,  and 
*b  functions  (usually  polynoaials)  of  x  and  y.  The  assumption  of 
linear  z  variation  of  inplane  strains  is  consistent  with  the  displace¬ 
ment  distribution  of  equations  (2.2).  However,  note  that  the  strains  in 
equations  (2.10)  are  not  directly  related  to  the  displaceaent  deriva¬ 
tives  (i.e.  not  related  to  £p*  in  equation  (2.3)). 

The  inplane  stresses  in  ltytf  i  ctn  bt  nltttd  to  ^  vit  substitute 
ing  squat ion  (2.10)  into  ths  inplant  strsss-strain  (plans  strsss)  rsla- 
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(2.11) 


The  matrix  Q1  ia  tha  inverae  of  tha  inplana  portion  of  /S1  aa  will  ba 
diaeuaaad  in  aora  datail  latar.  Nota  that  tha  inplana  atraaaaa  in  aqua¬ 
tion  (2.11)  vary  linaarly  in  z.  Tha  tranavaraa  ahaar  atraaaaa  and  tran- 
avaraa  normal  atraaa  in  layar  i  can  ba  ralatad  to  £  by  integrating  tha 
homogenaoua  equilibrium  aquationa;  i.a. 


(2.12a) 


(2.12b) 


(2.12c) 


In  general,  atraaaaa  from  aquation  (2.11)  are  aubatitutad  into  aquationa 
(2.12).  Tha  conatanta  of  integration  for  each  layar  in  aquationa  (2.12) 
are  uaad  to  aatiafy  traetion-fraa  conditiona  on  the  lower  aurfaca  (con- 
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stants  for  layer  1),  and  continuity  of  axz,  0yZ,  and  az  at  inter¬ 
layer  surfaces  (i.e.  at  z  m  h£  using  the  constants  of  layer  i  for  i  • 
2,3... I).  Details  of  these  calculations  are  involved  and  will  be  pre¬ 
sented  later  in  conjunction  with  an  8-node  element. 

Using  equations  (2.11)  and  (2.12),  and  the  general  procedure  just 
outlined,  it  will  be  possible  to  relate  the  stresses  in  layer  i  to  8  in 
the  fora 


(2.13) 


Note  that  the  stresses  satisfy  equilibrium  as  required  in  a  hybrid- 

stress  element. 

Equations  (2.9),  (2.13),  and  the  interpolation  for  wo  are  substi¬ 

tuted  into  the  hybrid-stress  functional  of  equation  (2.1),  and  the  fol¬ 
lowing  laminate  matrices  are  defined 


H  =  Z  Hi 

i 


(2.14a) 


& 


(2.14b) 


where,  invoking  the  mapping  of  equation  (2.5)  to  the  system 


y*  -  f  (  £  ££ 

-l  Cl 


(2.15a) 


-  /'/’  tT  51 


(2.15b) 


£  1  (  (  H*  f  lJl  J5 


(2.15c) 


In  aquation*  (2.15),  ii  ia  assumed  that  the  layer  stress  interpolation 
matrix,  P*,  and  the  prescribed  distributed  load,  p,  are  expressed  as 
functions  of  £  and  v  using  equations  (2.5).  Also,  |j|  is  the  Jacobian 
of  the  coordinate  transformation.  Equation  (2.1)  then  becomes 


ft  ATb 
■f  / 


ft.  -ft  &  $. 


(2.16) 


where  F  is  the  element  nodal  force  vector.  Equation  (2.16)  is  in  the 
standard  form  for  hybrid-stress  elements. 

He  also  consider  the  possiblity  that  additional  interrelations 
among  the  3  may  exist  which  can  be  expressed  in  the  general  form 
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-  O 


(2.17) 


Inter-relations  of  this  form  were  first  used  for  hybrid-stress  elements 
by  Mau  and  Pian  [6]  to  enforce  interlayer  traction  continuities  in  mul¬ 
tilayer  elements.  'hey  have  been  used  more  recently  as  an  alternate 
scheme  to  enforce  equilibrium  in  hybrid-stress  elements  [19-21] ,  and  to 
reduce  stress  fields  in  single  layer  hybrid-stress  Mindlin-type  plate 
elements  [18] .  They  could  also  be  used  to  enforce,  for  example,  partial 
equilibrium.  In  the  present  multilayer  elements,  the  inter-relations 
among  3,  expressed  by  equation  (2.17),  will  be  used  to  reduce  the  stress 
field  to  a  desirable  form  and  to  satisfy  the  zero  transverse  shear 
stress  conditions  at  the  upper  laminate  surface. 

Following  Reference  [6],  equation  (2.17)  is  introduced  into  equa¬ 
tion  (2.16)  via  the  Lagrange  Multiplier  method.  From  the  stationary 
condition  of  the  functional  it  can  be  shown  that  the  stress  parameters 
are  related  to  nodal  degrees  of  freedom  by 


A  -  a"  6* 


(2.18) 


and  the  element  stiffness  matrix  is  given  by 


k  -  gTi/"S 


(2.19) 


A.1  , 

where  H  1  is  the  augmented  matrix 
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( H',Sry>  K  H 


(2.20) 


The  formulation  just  presented  produces  a  multilayer  plate  element 
having  nodal  degrees  of  freedom  and  stress  parameters  which  are,  in  num¬ 
ber  ,  independent  of  the  number  of  layers.  Per  the  studies  of  References 
[9,  10],  a  multilayer  element  of  this  type  should  provide  reasonable 

predictions  for  thin  to  moderately  thick  laminates;  i.e.  for  typical 
thickness  to  span  ratios  of  h/L  “0.1  and  lower. 

In  early  phases  of  the  present  effort,  8~node  elements  of  this  type 
were  formulated.  Despite  the  fixed  number  of  stress  parameters  and  no¬ 
dal  degrees  of  freedom,  it  was  found  that  computation  times  to  form  k 
grew  excessively  with  increasing  numbers  of  layers.  This  was  found  to 
be  due,  to  a  large  extent,  to  the  complexity  of  ax2l ,  <JyZ1' •  211(1 

crzl  from  equations  (2.12)  and  their  contributions  to  H,  and  G  (equa¬ 
tions  (2.15)).  As  a  result,  it  was  decided  to  restrict  attention  to 
thin  multilyer  laminates  (i.e.  h/L  i  .05),  for  which  most  applications 
were  intended,  and  to  pursue  a  reduced  formulation  similar  to  that  re¬ 
cently  used  for  single  layer  plates  by  Spilker  [18]  which  allowed  for 
significant  savings  in  computation  times  for  k.  The  reduction  scheme 


and  corresponding  element  matrices  are  described  in  the  next  subsection. 
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2.3  STIFFNESS  MATRIX  FOR  A  THIN  MULTILAYER  PLATE  ELEMENT 

In  a  recent  study  [18] ,  Spilker  has  shown  that  a  hybrid-stress  for¬ 
mulation  for  single  layer  Mindlin-type  isoparametric  plate  bending  ele¬ 
ments  can  be  reduced  to  provide  computationally  efficient  isoparametric 
thin  plate  elements.  The  elements  utilize  independent  interpolations 
for  transverse  displacement  and  cross-section  rotations,  and  include  all 
components  of  stress.  It  was  shown  that  the  contributions  of  transverse 
shear  stress  and  transverse  normal  stress  to  the  internal  complementary 
energy  (i.e.  the  £  matrix)  could  be  neglected  for  thin  plates.  The  re¬ 
sulting  elements  were  found  to  converge  to  classical  thin  plate  theory 
solutions.  This  reduction  also  allowed  use  of  an  alternate  algorithm 
for  stiffness  matrix  calculation  which  produced  significant  savings  in 
computation  time.  This  reduction  is  developed  for  thin  multilayer 
plates  in  this  section,  and  an  efficient  implementation  for  stiffness 
matrix  calculation  is  given. 

2.3.1  Definition  of  the  Reduced  Functional 

In  order  to  define  a  reduced  hybrid-stress  functional  for  thin  mul¬ 
tilayer  plates,  it  is  necessary  to  consider  the  functional  of  equation 
(2.1)  along  with  the  assumptions/  interpolations  of  the  previous  sec¬ 
tion.  Of  particular  interest  are  the  through-thickness  distributions  of 
stresses  v1  and  strains  from  displacements,  c,1,  and  contributions  to 
each  term  in  equation  (2.1)  which  might  be  neglected  for  thin  plates.  A 
functional  will  then  be  defined  which  is  valid  for  thin  multilayer 
plates  and  a  computationally  efficient  implementation  of  this  functional 
will  be  presented. 
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Recall  that  inplana  strains  £p*  of  aquation  (2.3)  are  of  order 
z,  that  the  transverse  shear  strains  of  equation  (2.4)  are  con¬ 
stant  in  z,  and  that  the  transverse  normal  strain  £zl  "  0  for  all 
layers.  Per  equations  (2.10)  and  (2.11),  the  inplane  stresses  op* 
(i.e.  op*T*[ox*  ay*  ffxy^  )  are  °*  0Td*r  z«  Then  via 
equations  (2.12),  the  transverse  shear  stresses  (i.e.  j?t^  “ 
[crxz*  cr^1])  will  be  of  order  z^,  and  the  transverse  normal 
stress,  <rzS  will  be  of  order  z^. 

Now  consider  each  term  in  equation  (2.1)  as  plate  thickness  is  made 
small.  The  final  term  corresponds  to  nodal  loads  due  to  distributed  ap¬ 
plied  loads  and  is  not  affected  by  plate  thickness.  Consider  the  second 
term  in  equation  (2.1).  Contributions  of  the  product  £p*Tip*» 
integrated  through  the  thickness  will  be  of  order  h^  (for  present  pur¬ 
poses,  h  may  be  considered  as  a  normalized  thickness;  i.e.  plate  thick¬ 
ness  divided  by  a  typical  spanwise  dimension).  Likewise  the  contribu¬ 
tions  to  the  second  term  in  equation  (2.1)  of  the  product 
(transverse  shea-  stresses/atrains)  will  be  of  order 
h^.  The  product  oz*T?z*  does  not  contribute  since  «z*-0. 

In  view  of  the  above  contributions,  no  reduction  of  the  second  integral 
in  equation  (2.1)  is  possible  when  h  is  small. 

Reductions  of  the  first  integral  of  equation  (2.1)  are  possible  for 
small  h.  For  a  layer  having  a  general  fiber  orientation,  the  material 
property  matrix  will  have  the  following  form: 
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(2.21) 


^3} 

Expressions  for  the  are  given,  for  example,  in  Reference  [10]  . 

To  examine  the  effects  of  h,  it  is  convenient  to  expand  the  first  term 
in  equation  (2.1)  as  follows: 
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where 
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Referring  to  the  orders  of  h  in  equation  (2.22),  the  first  tern  is 
of  order  h^,  the  second  and  third  terns  are  of  order  h^,  the  fourth 
tern  is  of  order  h^,  and  the  fifth  tern  is  of  order  h^.  In  addi¬ 
tion,  each  of  the  second  through  fifth  terns  is  related  to  the  sane  set 
(or  a  subset)  of  3  as  the  first  tern.  Thus,  for  snail  h  the  second 
through  fifth  terns  in  equation  (2.22)  can  be  neglected  conpared  with 
the  first  tern. 

It  should  be  noted  that  it  is  possible  to  further  divide  the  pre- 
ceeding  discussion  into  contributions  due  to  stretching  and  those  due  to 
bending.  Although  the  contributions  of  stretching  and  bending  differ  by 
a  factor  h^,  the  relative  contributions  fron  each  of  the  terns  in 
equation  (2.22)  is  as  described  above  and  the  conclusions  drawn  above 
are  valid. 

With  the  reductions  described  above,  the  hybrid  functional  of  equa¬ 
tion  (2.1)  can  be  rewritten  for  thin  aultilayer  plates  in  the  following 
fora: 
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(2.24) 
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In  fuaaary,  the  functional  of  aquation  (2.24)  for  thin  multilayer  lami¬ 
nates  differs  from  the  general  laminated  plate  functional  in  equation 
(2.1)  in  that  contributions  of  transverse  shear  stresses  and  transverse 
normal  stress  to  the  internal  complementary  energy  have  been  neglected. 
This  is  analogous  to  the  reduction  performed  for  single  layer  plate  ele¬ 
ments  in  Reference  (18}.  Of  significance  in  practice  is  the  fact  that 
it  is  exactly  those  eliminated  terms  which  contributed  significantly  to 
stiffness  computation  time  in  the  non-reduced  element.  As  a  result  of 
this  reduction,  a  modified,  more  efficient  algorithm  for  stiffness  cal¬ 
culation  can  be  used.  This  algorithm  is  described  next. 

2.3.2  Definition  of  the  Element  Stiffness 

An  efficient  algorithm  can  be  defined  for  stiffness  formation  for  a 
thin  multilayer  plate  element  baaed  on  the  reduced  functional  in  equa¬ 
tion  (2.24).  He  begin  by  assuming  that  the  same  order  polynomial  is 
used  to  interpolate  each  component  of  inplane  strain,  both  stretching 
and  bending  contributions (see  equation  (2.10));  i.e. 

C*  =  £  A  +  }  &  A  *■ 

3  ^  (2.25) 

€*?  '  ~  /i5  r  h  ~  far 

where  a  is  a  vector  of  polynomials,  e.g. 


Using  the  notation  of  equation  (2.10) 
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(2.27b) 

The  inplane  stresses  in  layer  i,  £pl,  can  be  related  to  through 
equation  (2.11).  Note  that  siodifi cat ions  to  the  strain  distribution  in 
equations  (2.25)  can  be  achieved  by  introducing  the  necessary  interrela¬ 
tions  in  the  fora  of  equation  (2.17). 

The  transverse  shear  stresses  do  not  contribute  to  the  first  inte- 
gral  in  equation  (2.24)  (which  will  lead  to  the  H  matrix) ,  but  do  con¬ 
tribute  to  the  second  integral  in  equation  (2.24)  (which  will  lead  to 
the  G^aatrix).  To  define  £jl  for  layer  i,  equations  (2.25)  are  sub¬ 
stituted  into  equation  (2.11),  and  the  result  substituted  into  equations 
(2.12a)  and  (2.12b).  The  constants  of  integration  for  layer  1  are  used 
to  enforce  the  zero  traction  condition  at  the  lower  surface;  the  con¬ 
stants  of  integration  for  layer  i  are  used  to  satisfy  the  traction  con- 
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tinuity  at  z  •  hj  between  layer*  i  and  i-1  (i>l) .  After  sone  Manipu¬ 
lation,  it  can  be  ehown  that  <rxzl  and  are  given  by 
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where  the  coefficient  vector*  are  function*  of  x,y,  and  z,  and  are  given 
by 
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Inplane  strease*  £p*  and  transverse  shear  stresses  o^  are 

substituted  into  the  functional  of  equation  (2.24).  Also,  and 


<_(• 1  can  be  related  to  nodal  degree*  of  freedoa,  as  described  ear¬ 
lier,  and  substituted  into  equation  (2.24).  With  these  substi tutions , 
of  equation  (2.24)  can  be  put  in  the  fora 


■“-* =  %  "£rn  +  } 


(2.30) 


as  in  equation  (2.16). 

For  the  present  thin  aultialyer  plate  eleaent,  however,  the  utri- 
ees  H  and  G  can  be  defined  in  a  aore  convenient  fashion.  To  facilitate 
these  definitions,  we  assuae  J3  are  ordered  as  in  equation  (2.27b)  and  q, 
are  ordered  as  in  equation  (2.7).  We  then  define  the  following  inte¬ 
grated  vector  products  (note  that  the  integrals  are  perforaed  in  £,r;  and 
dA- | J | dfdv) : 


a?  a. 


(2.31) 
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In  addition,  we  define  the  following  "material  parameter*"  which  come 
about  in  the  through~thickness  integrations,  summed  over  all  layers: 


r 

t 
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(2.33a) 


(2.33b) 
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With  these  definitions,  H  is  given  by 
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where  E^j  are  terms  in  the  symmetric,  positive  definite  matrix  E  given 
by 
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The  inverse  of  H  (needed  for  k)  is  then  given  by  [23] : 
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The  matrix  G  can  be  defined  as 
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Note  that  the  ordering  of  (equation  (2.7))  used  to  define  £  is  conven¬ 
ient  but  not  conventional.  In  practice,  q  will  be  ordered 

Ltto,  -rfe/  flot  &*!  &yi  y«2  ^az.  ^4  ■'•  ^ 

$  (2.38) 


and  terns  in  C  fron  equation  (2.37)  will  be  assigned  to  the  proper  col¬ 
umn  location  corresponding  to  equation  (2.38). 

In  most  applications  of  the  present  thin  multilayer  plate  formula¬ 
tion,  additional  interrelations  among  the  in  the  form  of  equation 
(2.17)  will  be  needed.  For  example,  the  stress  field  resulting  from  the 
strain  field  of  equations  (2.25)  will  contain  complete  polynomials  in 
order  to  preserve  natural  invariance  [24] .  The  resulting  thin  plate 
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constraints  (to  be  discussed  in  the  next  section)  may  be  excessive,  and 
thus  a  reduction  of  the  stress  field  (which  preserves  invariance)  may  be 
necessary.  Such  a  reduction  could  be  put  in  the  form  of  equation 

(2.17).  As  a  second  example,  it  is  noted  that  the  stress  field  used  for 
the  present  thin  multilayer  plate  elements  satisfies  equilibrium,  satis- 
fies  lower  laminate  surface  traction-free  conditions,  and  satisfies 
traction  continuity  on  interlayer  surfaces.  Recall  that  the  preceeding 
traction  conditions  were  satisfied  using  the  constants  of  integration  in 
equations  (2.12a)  and  (2.12b).  If  the  inplane  tractions  on  the  upper 
laminate  surface  (i.e.  <7XZ  and  <7yZ)  are  to  be  set  to  zero,  this  up¬ 
per  surface  condition  must  be  enforced  by  introducing  interrelations 
among  the  in  the  form  of  equation  (2.17).  These  two  kinds  of  stress 
field  reductions  will  be  illustrated  later  in  the  development  of  an 
8-node  element. 

With  H,  G,  and  jl,  defined  as  just  described,  the  element  stiffness 
can  be  defined  by  equation  (2.19).  The  major  computational  advantages 
brought  about  by  the  reduced  formulation  are  summarized  in  the  follow¬ 
ing: 

(a)  Perhaps  the  most  significant  computational  saving  results  because 
the  only  layer  dependent  operations  in  the  formation  of  are  the  defi¬ 
nitions  of  the  constants  in  equations  (2.33).  These  operations  are  so 
insignificant  that  there  is  no  real  increase  in  computation  time  from  1 
to  100  layers.  In  effect,  computation  time  is  independent  of  the  number 


of  layers. 
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(b)  The  area  integrals  in  equations  (2.31)  and  (2.32)  can  be  efficiently 
evaluated  numerically.  Only  the  lower  triangle  of  (equations  (2.31)) 
need  be  evaluated.  Also,  if  desired,  the  integrals  could  be  evaluated 
analytically  since  no  |j|  will  occcur  in  the  denominator. 

(c)  The  calculation  of  would  normally  involve  a  system  of  the  or¬ 
der  of  the  number  of  8.  Instead,  the  definition  of  IJ  '  by  equation 
(2.36)  requir*  inversion  of  a  symmetric  6x6  matrix  (E)  and  inversion  of 
a  symmetric  matrix  (d>)  which  involves  approximately  one-sixth  the  number 
of  6  s. 

(d)  The  "cost"  of  items  (b)  and  (c)  is  the  more  complicated  stiffness 
expression  in  equations  (2.19)  and  (2.20).  Since  the  constraints  in 
equation  (2.17)  do  not  involve  equilibrium,  they  are  generally  few  in 
number  compared  with  the  number  of  8's.  Thus  £  contains  many  zeroes, 
and  operations  involving  It  can  be  streamlined  to  avoid  most  multiplica¬ 
tions  by  zero. 

Illustrative  comparisons  of  stiffness  computation  times  for  an 
8-node  element  will  be  given  later. 

2.3.3  Constraints  in  the  Thin  Plate  Limit 

In  the  earlier  study  of  single  layer  hybrid-stress  plate  elements 
[19] ,  it  was  shown  that  the  reduced  thin  plate  element  was  identical  to 
the  nonreduced  (moderately  thick)  plate  element  in  the  thin  plate  limit. 
Thus,  the  performance  of  the  thin  plate  reduced  element  will  be  guided 
by  the  thin-plate~limit  constraints  of  the  nonreduced  element. 


The  sane  holds  for  the  present  multilayer  elenents:  the  reduced 

thin  multilayer  plate  element  will  be  identical  to  the  nonreduced  multi¬ 
layer  element  in  the  thin  plate  limit.  The  thin-plate-limit  constraints 
on  the  nonreduced  element  are  enforced  for  all  plate  thickneses  on  the 
reduced  thin  plate  element.  If  we  denote  by  jcTt  the  transverse  shear 
stresses  integrated  through  the  thickness  and  summed  over  all  layers, 
then  the  thin-plate-limit  constraints  of  the  nonreduced  formulation  are 


^  -p 

f  £*  JA  »  O  Jn  ££  jo  (2.39) 

N J  ^ 


Equation  (2.39)  represents  weighted,  integrated,  Kirchhoff-type  con¬ 
straints.  The  form  is  analogous  to  the  constraints  for  the  single  lay¬ 
ers  elements  of  References  [15-18] .  Thus,  by  analogy  with  the  single 
layer  elements,  multilayer  elements  of  the  present  type  which  are  based 
on  a  fixed  set  of  (laminate)  stress  parameters  should  have  spanwise 
distributions  of  a^z  and  which  are  identical  to  those  of  the 

corresponding  non-locking  single  layer  hybrid-stress  plate  elements 
[19] . 

In  the  reduced  multilayer  thin  plate  element,  the  constraints  of 
equation  (2.39)  are  imposed  for  all  plate  thicknesses.  If  the  nonre¬ 
duced  element  is  nonlocking  in  the  thin  plate  limit,  the  reduced  element 
will  converge  to  the  lamination  (thin  plate)  theory  solution  regardless 
of  plate  thickness. 
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2.4  FORMATION  AND  EVALUATION  OF  AN  INVARIANT  8-NODE  THIN  MULTILAYER 
PLATE  ELEMENT 

2.4.1  Element  Displacement  and  Stress  Interpolations 

To  illustrate  the  formulation  of  a  thin  multilayer  plate  element 
based  on  the  reduced  formulation  of  the  proceeding  section,  we  consider 
an  8-node  multilayer  plate  element  as  shown  in  Figure  2.  The  8-node  el¬ 
ement  has  been  selected  in  view  of  studies  of  single  layer  hybrid-stess 
(Hindi in-type)  plate  elements,  wherein  it  is  shown  that  a  nonlocking, 
naturally  invariant  8-node  element  of  correct  rank  can  be  developed 
[18].  That  element,  denoted  QH3  in  Reference  [18],  has  served  as  the 
guide  for  a  multilayer  element  with  independent  layer  displacement  and 

I 

stress  fields  [11]  and  a  reduced  thin  single  layer  element  [19] ,  and 
will  guide  the  development  of  the  present  8-node  thin  multilayer  plate 

element. 

The  displacement/rotation  interpolation,  and  the  geometric  mapping 
(since  the  element  is  isoparametric)  used  for  the  8-node  element  is  the 
standard  8-node  biquadratic  Serendipity  shape  functions  (expressed  in 
the  system) .  The  shape  functions  N£(|,n),  i  “  1,2... 8  for  equa¬ 
tions  (2.5)  through  (2.7)  can  be  found  in  most  texts  on  finite  element 
methods.  With  five  degrees  of  freedom  per  node,  the  present  element  has 
40  dof,  independent  of  the  number  of  layers. 

The  stress  field  for  our  model  single  layer  plate  element,  QH3,  was 
obtained  by  starting  with  complete  cubic  polynomials  in  x,y  for  each  of 
the  inplane  stresses.  The  stress  field  was  then  subjected  to  the  con¬ 
straint  that  ?z  be  independent  of  x,y. 

To  produce  the  analog  multilayer  thin  plate  element,  the  strain  in¬ 
terpolation  vector,  a,  in  equations  (2.25)  should  be  a  complete  cubic 
polynomial,  i.e. 


As  a  rssult,  tha  unconstrained  stress  £ield  contains  a  total  of  60  3  ' s, 
with 
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Given  the  N;(£,tj)  and  a,  the  H  and  G  matrices  can  be  formed. 

'Vj  ** 

The  60  @  stress  field  resulting  from  substituting  equation  (2.40)  into 
equation  (2.25)  must  be  reduced  to  correspond  to  QH3.  In  particular, 
the  distribution  of  ~cz  given  by 
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**  stt 

must  be  independent  of  x  and  y.  Since  ?z  is  linear  in  x  and  y,  this 
produces  two  constraints  as  follows: 
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In  addition,  we  require  that  <?xz  and  be  zero  at  the  upper 

laminate  surface.  Since  ox z  and  & are  complete  quadratic  polyno- 
mials,  this  produces  a  total  of  12  constraints,  as  follows: 
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Equations  (2.43)  and  (2.44)  represent  14  interrelations  among  the 
3’s  which  can  be  put  in  the  £ora  of  equation  (2.17)  and  thus  define  the 
R  matrix.  If  we  consider  these  interrelations  as  reducing  the  effective 
number  of  3's,  we  can  say  that  the  present  element  has  a  total  of  46 
3's.  In  future  discussions,  we  will  refer  to  this  8-node  thin  multilay¬ 
er  plate  element  as  MQH3T. 

Element  MQH3T  possesses  correct  rank,  independent  of  the  number  of 
layers.  It  is  also  naturally  invariant.  As  it  is  based  on  the  reduced 
formulation  of  the  previous  section,  it  should  be  expected  that  MQH3T 
will  converge  to  lamination  theory  solutions,  and  that  convergence  be¬ 
havior  will  be  like  the  single  layer  element  QH3  [18] . 

2.4.2  Element  Stiffness  Computation  Time 

Computation  time  (on  an  IBM  370/168  using  double  precision)  for  the 
stiffness  matrix  of  MQH3T  is  3.38  seconds,  essentially  independent  of 
the  number  of  layers.  A  breakdown  of  that  time  is  given  in  Table  1. 
The  largest  single  block  of  time  is  for  the  final  matrix  multiplica¬ 
tions.  The  next  largest  block  is  the  time  to  compute  H“*. 


Wv  V  *>  V  ■ 


TABLE  1 


Computation  Time  for  the  Stiffness  Matrix  of  Element  MQH3T 


OPERATION 

CPU  (secs) 

Zero/Initialize  Matrices 

0.034 

Evaluate  area  Integrals  in  eqns. 

(2.31)  and  (2.32)  (Needed  for  Hh 
and  G)  Using  a  4x4  Gauss  Rule 

0.124 

Form  <tr^  (eqn.  (2.31)  and  E-* 

(eqn. (2.35)) 

0.012 

Form  H~1  (eqn. (2.20)) 

0.843 

Form  G  (eqn.  (2.37)) 

0.559 

Compute  k“G^H~*G 

1.810 

TOTAL 

3.382 

Some  comparisons  with  computation  times  for  other  element  stiff¬ 
nesses  are  useful.  Early  attempts  to  produce  a  non-reduced  8-node  mul¬ 
tilayer  element,  hereafter  denoted  MQH3,  having  fixed  8's  and  degrees  of 
freedom  resulted  in  stiffness  computation  times  of  27.9  secs  for  1  lay¬ 
er,  and  109. A  secs  for  3  layers.  This  is  what  motivated  the  study  of  a 
reduced  formulation.  The  8-node  multilayer  element  MLPQH3  of  Reference 
[11]  having  independent  layer  stress  fields  and  independent  layer  dis¬ 
placement  fields  requires  30.4  secs  for  3  layers  (note — less  CPU  time 
than  MQH3) .  Finally,  the  single  layer  hybrid-stress  8-node  plate  ele¬ 
ment  (bending  only)  based  on  a  reduced  formulation,  QH3T,  of  Reference 
[18]  requires  1.8  secs.  These  comparisons  illustrate  the  computational 
advantages  of  MQH3T. 
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2.4.3  Examp 1 e  Problems  and  Numerical  Results 

2.4.3. 1  Cylindrical  Bending  of  a  3-layer  90/0/90  Cross-ply 
The  problem  considered  is  a  3  layer  90/0/90  cross-ply  ( 6  measured 
from  the  x-axis)  plate  of  infinite  length  in  the  x-direction,  length  L 
in  the  y-direction,  simply  supported  at  y*Q,L  and  subject  to  a  sinusoi¬ 
dal  distributed  transverse  load  (uniform  versus  x) .  Mesh  and  boundary 
conditions  are  given  in  Figure  3.  Each  layer  is  of  equal  thickness, 
h/3,  and  layer  material  properties  are  E11  ■25x10*  psi.,  E22“10‘P*i*» 
*'^2“*25,  «'23“*25,  Gj2“»5xl0*  psi.,  G23“.2xl0*  psi.  This  problem 
has  been  selected  because  an  elasticity  solution  exists  which  is  valid 
for  all  laminate  thickness  ratios  [2^ . 

Quantities  of  interest  are  the  transverse  displacement  w,  the 
through  thickness  distribution  of  in-plane  normal  stress  oy  at  the 
plate  center  (y-L/2) ,  and  transverse  shear  stress,  o^,  at  the  bounda¬ 
ry  (y-0) . 

Results  will  be  presented  in  terms  of  the  following  normalized 
quantities . 
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The  quantity  qg  is  the  magnitude  of  the  distributed  sinusoidal  load. 
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Figure  3  P rob lea  definition  and  nash  for  tha  cylindri 
banding  of  a  seei-inf inite  plata  under  sinusoidal  loading 
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For  a  thin  laminate  the  linear  through  thickness  inplane  displace* 
aent  field  assuaed  in  MQH3T  is  adequate.  Convergence  of  the  center 
transverse  displacement  normalized  by  the  exact  value,  [24]  for  a  thin 
(h/L“0.01)  laminate  is  shown  in  Figure  4.  The  behavior  of  the  multilay¬ 
er  composite  plate  element  with  complete  stress  contributions  (MQH3)  is 
compared  with  the  reduced  thin  aultilayer  plate  element  (MQH3T)  which 
excludes  the  contributions  of  transverse  shear  stresses  and  transverse 
noraal  stress  to  the  internal  coapleaentary  energy  (i.e.  the  £  matrix) . 
The  two  elements  converge  in  a  similar  manner  and  both  give  very  good 
results  with  element  MQH3  converging  only  slightly  more  quickly  than 
MQH3T:  this  is  because  the  results  are  normalized  by  the  exact  solution, 
rather  than  by  lamination  theory. 

For  thin  laminates,  the  performance  of  element  HQH3T  is  essentially 
identical  to  the  element  MQH3.  As  thickness  ratio,  h/L,  increases  the 
more  general  element  MQH3  is  necessary  to  produce  accurate  predictions. 
To  compare  the  range  of  applicability  of  MQH3T  and  MQH3,  plates  having 
thickness  ratios  in  the  range  h/L“0.01  to  h/L«0.25  have  been  analyzed 
using  a  three  element  aesh  (converged  per  Figure  2).  The  results  ob¬ 
tained  for  the  center  transverse  displaceaent ,  "vTc  are  given  in  Table  2 
for  element  MQH3  and  in  Table  3  for  eleaent  MQH3T.  Element  MQH3  produc¬ 
es  accurate  predictions  of  w  over  a  wide  range  of  thickness  ratios  from 
thin  (h/L**0.01)  to  thick  (h/L~0.25)  laminates.  In  Table  3,  the  predict¬ 
ed  «TC  using  element  NQH3T  does  not  increase  with  increasing  h/L.  This 
is  consistent  with  laaination  theory  (as  expected),  which  gives  poor 
predictions  of  center  transverse  displaceaent  for  values  of  h/L>.01. 
Note  that  laaination  theory  produces  a  value  of  "w". 50966  for  this  prob¬ 
lem,  in  essential  agreeaent  with  the  MQH3T  results. 


Number  of  elements 


Figure  A  Convergence  of  the  center  transverse  displacement 
normalized  by  the  elasticity  solution  for  the  cylindrical 
bending  of  a  plate  for  which  h/L-0.01. 


table  2 


Effects  of  thickness  ratio,  h/L,  on  the  predicted  center  transverse 
displacement  of  element  MQH3  for  the  cylindrical  bending  of  a  90/0/90 

cross-ply  for  a  3  element  mesh 


h/L  ( 

Mc> exact 

(“^predicted 

predicted/exact 

0.25 

3.0233 

2.8872 

1.047 

0.10 

0.9278 

0.9333 

0.994 

0.05 

0.6160 

0.6167 

0.999 

0.0375 

0.5689 

0.5702 

0.998 

0.01 

0.50906 

0.5133 

0.992 

TABLE  3 

Effects 

displacement 

of  thickness  ratio,  h/L,  on  the  predicted  transverse 
of  element  MQH3T  for  the  cylindrical  bending  of  a  90/0/90 
cross-ply  using  a  3  element  mesh 

h/L 

(wc5 exact 

(wc) predicted 

predicted/exact 

0.25 

3.0233 

0.50104 

0.1657 

0.10 

0.9278 

0.50123 

0.5402 

0.05 

0.616 

0.50123 

0.8137 

0.0375 

0.5689 

0.50122 

0.8811 

0.01 

0.50906 

0.50123 

0.9846 

In  the  range  h/L>.01,  MQH3  is  clearly  the  sore  accurate  element 
since  the  contribution  of  transverse  shear  stresses  and  transverse  nor¬ 
mal  stresses  to  the  internal  complementary  energy  have  been  retained. 
However,  as  illustrated  earlier,  computation  time  for  k  of  MQH3  grows 
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excessively  with  increasing  numbers  of  layers.  Thus,  use  of  MQH3  for 
thin  plates  is  not  appropriate.  The  reduced  element  is  clearly  the  bet¬ 
ter  element  for  h/L<.01  because  of  the  significant  savings  in  CPU  time 
and  the  fact  that  there  is  no  loss  of  accuracy  within  the  range  of  h/L 
to  be  considered  here  (thin  plates).  Therefore  only  the  reduced  element 
will  be  used  for  subsequent  tests. 

The  prediction  of  spanwise  distributions  of  stresses  can  be  illus¬ 
trated  using  the  90/0/90  example  problem.  Spilker  et  al.  [23]  estab¬ 
lished  the  existence  of  optimal  sampling  points  for  hybrid-stress  ele¬ 
ments.  Briefly,  approximate  strains  approach  the  exact  strains  £u, 
(where  u  are  the  exact  displacements  and  D  is  the  differential  operator 
matrix  corresponding  to  the  linear  strain-displacement  equations)  in  a 
weighted  least  squares  sense.  Owing  to  a  property  of  Gauss  quadrature, 
the  Gauss  stations  corresponding  to  the  Gauss  rule  which  just  integrates 
the  stress  field  are  the  optimal  sampling  points. 

To  examine  the  predicted  spanwise  distributions,  it  is  convenient 
to  consider  the  (normalized)  moment  (from  cy)  and  shear  force  (from 
axz)  distributions  given  by: 

h 

*  (2.46) 
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The  present  element  MQH3T  is  based  on  a  single  layer  element  QH3 
[17] .  Evaluation  of  the  single  layer  elenrnt  determined  that  conver¬ 
gence  of  M  at  nodes  is  in  general  somewhat  better  than  at  the  2x2  Gauss 
stations  while  the  results  obtained  for  Q  at  the  2x2  Gauss  stations  are 
extremely  accurate.  Using  element  QH3  as  a  guide,  results  are  plotted 
in  Figures  5  and  6  for  M  and  ‘Q  for  a  thin  90/0/90  plate,  h/L“.01  and  us¬ 
ing  3  elements.  Results  for  M  are  plotted  versus  y  along  the  left  edge 
of  the  mesh.  The  prediction  of  TT  is  reasonably  accurate  at  the  nodes 
and  center  of  each  element.  Results  for  Q  are  plotted  versus  y  along  a 
line  which  passes  through  the  2x2  Gauss  stations.  The  predicted  distri¬ 
bution  of  Q  is  poor  on  a  pointwise  basis  but  the  optimal  sampling  points 
corresponding  to  the  2  point  Gauss  rule  are  evident. 

The  through  thickness  distributions  of  Oy  at  y~L/2  and  transverse 
shear  stress  near  y-0  for  thin  laminated  (h/L-0.01)  are  shown  in 
Figures  (7a~7b) .  A  3  element  mesh  is  used  and  the  lamination  theory  so¬ 
lution  is  also  shown.  The  predicted  distribution  of  <7y  at  y/L*0.5 
agrees  exactly  with  the  lamination  theory  solution  (Figure  7a).  The 
predicted  distribution  of  ayz  calculated  at  the  2x2  Gauss  station 
nearest  to  y“0  agrees  very  closely  with  lamination  theory.  Figure  7b. 

2. 4. 3.2  Angle  ply  Laminates  under  Transverse  Loading 

The  final  example  problem  considered  is  a  2  layer  laminate  of  angle 
ply  (10)  contruction,  subject  to  a  uniformly  distributed  transverse  load 
(see  Figure  8).  For  present  purposes,  only  square  plates  are  considered 
so  that  a“b  and  a  uniform  NxN  mesh  is  used;  note  that  the  entire  plate 
must  be  modeled  since  fiber-orientation  induced  stretehing/bending 


Figure  5  Normalized  moment  N  distribution  vs  y  for  a  thin 
(h/L"0.01)  90/0/90  cross-ply  in  cylindrical  bending  using  a 
3  element  mesh. 


°yz  *t  2x2  Gauss  stations  naar  y*0 


Fl*ur*  Through  thieknass  strass  distribution  for  thin 
(h/l-0.01)  90/0/90  eross-piy  laainatas  in  cylindrical 
banding  using  a  3  alaaant  nash. 


coupling  exists  which  eliminates  the  usual  (i.e.  single  layer)  x-y  plane 
symmetries.  Results  will  be  compared  with  a  series  solution  (20  terms 
used  here)  given  by  Whitney  [25]  since  that  solution  is  based  on  lamina¬ 
tion  theory. 

All  results  are  obtained  £or  a*b“10,  h**0.02,  so  that  h/L*.002. 

Layer  material  properties  are  Eu**40xl0*  psi.,  E22"^®‘  P»i»» 

vl2mv23m'^’  ^12"®23"'^xl0*  psi. 

For  all  cases,  a  uniformly  distributed  transverse  load  is  applied. 
Plates  are  considered  for  which  all  sides  are  simply  supported  and  for 
which  all  sides  are  ideally  clamped.  Quantities  of  interest  to  assess 
element  performance  are  transverse  displacement,  wc  and  moments,  Mx 
and  My  at  the  center  of  the  plate.  As  well  as  inplane  displacements  u 
and  v  along  the  edge  of  the  plate. 

Results  are  presented  first  for  simply  supported  plates.  Limited 
convergence  results  for  the  case  9m  15°  are  given  in  Table  4.  Table  4 
shows  all  quantities  of  interest  to  be  within  2.5  per  cent  of  the  exact 
value  for  a  6x6  mesh  of  the  entire  plate. 

Using  the  6x6  mesh  of  the  whole  square  plate,  the  effect  of  9  on 
the  prediction  of  w,  is  shown  in  Table  5.  Typically,  9  has  little  ef¬ 
fect  on  the  accuracy  of  the  predictions. 

The  second  example  problem  is  a  plate  clamped  on  all  sides.  Con¬ 
vergence  results  for  9  m  15°,  25s,  and  45°  are  given  in  Table  6.  The 
results  are  poor  for  small  9,  i.e.  9m  15°  and  good  for  9m  45°.  Table  6 
shows  the  transverse  displacement  for  the  case  9m  45°  to  be  within  7  per 
cent  of  the  exact  value  for  a  6x6  mesh. 


(a)  geometry 


N  equal  divisions 


(b)  mesh  and  boundary  conditions  (top  view) 


Figure  8  Problem  definition  and  mesh  for  a  2  layer 
rectangular  plate  of  angle-ply  (i 6 )  contstruetion  under 
uniform  transverse  loading. 


TABLE  4 


Convergence  of  s  simply  supported  square  angle  ply  plate  under  uniform 

transverse  load  t  =rLi5  * 


Kesh  (whole  plate) 

2x2  4x4  6x6 


U0  at  (a/2,0) 
exact*. 00771 

pred. 
pred ./ex. 

0.009 

1.171 

0.0079 

1.0295 

0.0078 

1.013 

wq  at  (a/2,b/2) 
exact*8.927 

pred. 
pred. /ex. 

8.113 

.9088 

8.9883 

1.007 

8.9486 

1.0024 

Hx  at  (a/2,b/2) 
exact"10.42 

pred. 
pred. /ex. 

8.6127 

0.754 

11.8877 

1.041 

11.6890 

1.0236 

v0  at  (0,b/2) 
exact*0. 03281 

pred. 
pred. /ex. 

0.03356 

1.0228 

0.03279 

0.996 

0.03279 

0.9994 

My  at  (a/2,b/2) 
exact*1.234 

pred. 
pred. /ex. 

1.1193 

.907 

1.241 

1.006 

1.2332 

0.9994 

TABLE  5 

Results  for  a  thin 
construction 

siaply  supported  square  plate  of  angle  ply 
under  transverse  load  for  a  6x6  aesh. 

Fiber  Orientation 

wc 

wexact 

predicted/exact 

±5° 

595.5 

592.0 

1.006 

±15° 

894.9 

892.7 

1.002 

±25° 

983.3 

983.8 

1.000 

±35° 

944.3 

945.1 

0.999 

±45w 

914.4 

915.2 

0.999 

Using  the  6x6  aesh  of  the  whole  square  plate,  the  effect  of  6  on 
the  prediction  of  the  transverse  displaceaent  is  shown  in  Table  7.  The 


TABLE  6 


Convergence  of  the  transverse  displacement  at  the  center  of  a  clamped 
square  angle-ply  plate  under  uniform  transverse  load 


Fiber  2x2 

Orientation  we 

±15° 

±25° 

±45° 


pred./ex. 

2.5378  1.2006 
3.3457  1.1365 
3.6611  1.0130 


Mesh 

4x4 

v£  pred./ex. 

2.4057  1.138 

3.0768  1.045 

3.5004  0.970 


6x6 

wc  pred./ex. 

2.5109  1.1879 

3.2143  1.0920 

3.8569  1.0680 


problems  occur  when  the  length  of  an  element  is  much  greater  in  one  di¬ 
rection  than  the  other.  As  has  been  stated  earlier  the  element  we  are 
now  considering  is  analogous  to  a  single  layer  element  QH3  of  Reference 
[17].  In  the  single  layer  element  study  [17],  these  effects  were  most 
apparent  for  aspect  ratios  of  a/b"  2  and  3  in  the  case  of  clamped  bound¬ 
aries,  where  as  much  a  5  per  cent  was  found.  Considering  the  present 
problem  of  laminated  plates  when  theta  is  small  and  the  fibers  nearly 
line  up  with  the  x-axis,  the  element  is  much  stronger  in  the  x-direction 
than  in  the  y-direction.  This  may  be  causing  a  "material  aspect  ratio" 
effect  analogous  to  the  "geometric  aspect  ratio"  effect  observed  in  iso¬ 
tropic  single  layer  plate  elements -  further  study  to  the  phenomenon  is 

needed. 

2.4.4  Concludina  Remarks 

This  chapter  has  dealt  with  a  hybrid  formulation  of  an  isoparame¬ 
tric  element  for  the  analysis  of  thin  multilayer  laminated  composite 
plates.  The  principal  objective  was  the  development  and  evaluation  of 
an  element  in  which  the  displacement  behavior  is  characterized  by  lami- 
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nate  reference  plane  inplane  and  transverse  displacements  and  laminate 
non-normal  cross  section  rotations  and  in  which  the  components  of  stress 
are  related  to  a  set  of  laminate  stress  parameters;  i.e.  so  that  there 
is  no  increase  in  the  number  of  degrees  of  freedom  or  the  number  of 
stress  parameters  as  the  number  of  layers  increases.  For  thin  multilay¬ 
er  laminates,  a  reduced  formulation  was  defined  which  allows  for  a  sig¬ 
nificantly  more  efficient  stiffness-generating  algorithm.  The  results 
of  this  sudy  have  been  applied  to  the  development  of  element  MQH3T. 

In  summary,  element  MQH3T  gave  good  convergence  performance.  The 
element  retains  a  high  degree  of  accuracy  for  thin  plates  having  a 
thickness  ratio  less  than  h/L"0.01.  The  spanwise  distribution  of  M 

pointed  out  that  the  beat  points  to  calculate  ffy  are  at  the  sides  and 

center  of  each  element.  The  spanwise  distribution  of  Q  indicated  that 
the  bast  points  to  calcuated  Qyz  are  at  the  2x2  Gauss  stations.  The 
poor  results  displayed  for  small  9  in  the  case  of  an  angle-ply  plate 
which  has  been  clamped  on  all  edges  are  thought  to  be  due  to  a  type  of 
"material"  aspect  ratio  problem.  For  the  case  of  an  angle-ply  plate 

which  is  simply  supported  on  all  edges  the  results  are  very  good  and  9 

has  little  effect  on  the  accuracy  of  the  predictions. 


Chapter  III 

FORMULATION  OF  THE  CONTACT  PROBLEM 


3.1  INTRODUCTION 

The  finite  element  method  is  now  widely  recognized  as  a  powerful 
tool  in  solving  problems  in  structural  and  continuum  mechanics.  The  fi¬ 
nite  element  method  owes  its  popularity  largely  to  its  generality-  It 
permits  easy  modeling  of  arbitrary  shapes,  with  arbitrary  loads,  bounda¬ 
ry  conditions  and  can  include  the  effects  of  material  and/or  geometrical 
nonlinearities.  The  unique  characteristics  of  the  contact  problem  make 
an  analysis  difficult  even  by  the  finite  element  method.  As  the  load 
varies,  the  region  of  contact  between  bodies  changes  and  sliding  rela¬ 
tive  to  each  other  occurs.  This  sliding  may  be  accompanied  by  friction¬ 
al  forces  which  act  to  oppose  the  sliding  motion. 

The  work  done  thus  far  in  the  area  of  finite  element  analysis  of 
contact  problems  can  be  differentiated  in  the  following  areas.  First, 
the  formulation  is  based  on  either  the  displacement  model  or  on  various 
hybrid  and  mixed  models.  Second,  the  models  are  capable  of  considering 
contact  problems  in  which  there  may  be  node  to  element  side  contact  or 
they  are  capable  of  handling  only  problems  in  which  there  is  node  to 
node  contact.  The  last  distinguishing  characteristic  is  that  the  effect 
of  friction  along  the  contact  surface  may  or  may  not  be  included.  The 
largest  portion  of  the  work  done  has  utilized  the  displacement  model 
formulation  and  can  only  account  for  moderate  sliding  between  contacting 
bodies  (i.e.  node  to  node  contact). 


This  method  was  first  applied  by  Yamada,  et.al.  [26]  to  the  Hertz 
contact  problem  of  two  identical  cylinders  pressed  together.  The  meshes 
were  divided  so  that  nodes  of  one  cylinder  were  in  contact  with  the 
nodes  of  the  opposing  cylinder. 

White  and  Enderby  [27]  applied  the  finite  element  analysis  to  the 
problem  of  a  connecting  rod  eye  loaded  by  means  of  a  pin.  They  used 
linear  triangular  elements  to  model  the  two  deformable  bodies  and  spe¬ 
cial  overlapping  connecting  elements,  which  were  assigned  stiffness  to 
model  the  surface  of  contact.  Iterations  were  required  to  obtain  the 
correct  displacements . 

Gaertner  [28]  investigated  plane  elastic  contact  of  solids  subject 
to  small  strains  with  friction.  Linear  triangular  displacement  finite 
elements  were  used.  Node  to  node  contact  was  considered  and  an  itera¬ 
tive  method  was  utilized  to  find  the  equilibrium  state. 

Fredriksson  [29]  introduced  a  2-dimensional  contact  constitutive 
relation;  a  general  slip  criterion  with  associated  slip  rule  was  includ¬ 
ed.  The  contact  surface  was  limited  to  extending  in  the  node  to  node 
configuration.  In  order  to  reduce  the  size  of  the  problem,  degrees  of 
freedom  not  related  to  the  nonlinearities  were  eliminated  using  the  su¬ 
per  element  technique.  A  Coulomb  type  of  slip  criterion  with  an  associ¬ 
ated  slip  hardening  rule  was  assumed. 

Gap  elements  were  utilized  by  Stadter  and  Weiss  [30]  to  model  the 
gap  process.  The  solution  was  achieved  through  a  iterative  procedure 
which  adjusts  the  modulus  of  the  gap  elements.  The  procedure  was  ap¬ 
plied  to  the  thermostruetural  analysis  of  the  heat  shield  of  a  radioiso¬ 
tope  thermoelectric  generator  design. 
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The  method  presented  by  Okamoto  and  Nakazawa  [31]  dealt  with  the 
irreversibility  resulting  from  stick-slip  phenomenon.  The  loading  was 
incremental  where  the  magnitude  of  each  increment  of  load  was  calculated 
to  cause  a  change  in  one  node  pair  along  the  contact  surface.  The  con¬ 
tact  conditions  were  additional  conditions  independent  of  stiffness 
equations. 

Sachdeva  and  Ramakrishnan  [32]  generalized  a  procedure  developed  by 
Francaville  and  Zienkiewicz  [33]  to  include  frictional  effects  under 
proportional  loading.  The  method  uses  the  flexibility  matrix  which  was 
obtained  by  inversion  of  the  condensed  stiffness  matrix  in  which  all 
nodes  were  eliminated  <  'ept  those  where  contact  was  likely  to  take 
place. 

A  direct  automated  procedure  for  frictionless  contact  problems  is 
presented  by  Mahmoud  et.al.  [34].  This  method  circumvents  the  need  for 
the  inclusions  of  interface  elements  and  does  not  calculate  contact 
pressure.  Simplicity  seems  to  be  the  primary  goal  of  its  authors.  The 
method  applies  only  to  contact  problems  in  which  there  is  no  tangential 
motion  relative  to  the  contact  surface  and  to  contact  problems  which  are 
purely  advancing  or  receding.  The  direct  incremental  method  used  here 
means  the  applied  load  was  scaled  so  that  it  would  just  close  the  re¬ 
spective  gap  yet  not  induce  the  interpenetration  to  the  contact  surfac¬ 
es  . 

Displacement  models  have  also  been  developed  to  account  for  exten¬ 
sive  sliding  contact,  i.e.  nodes  of  one  body  may  come  in  contact  with 
the  side  of  the  element  on  the  boundary  of  the  opposing  body.  Therefore 
these  models  also  account  for  friction  along  the  contact  surface. 


Chan  and  Tuba  [35]  solved  the  problem  of  contact  using  linear  tri¬ 
angular  elements  derived  by  an  assumed  displacement  model.  The  model 
allowed  for  node  to  node  and  node  to  element  boundary  contact  and  cou¬ 
lomb  friction  betweed  the  two  bodies.  The  effects  of  clearance,  fric¬ 
tion,  and  load  on  stresses  in  turbine  blade  fastenings  were  studied  and 
compared  with  photoelastic  experimental  results. 

Large  displacement  contact-impact  problems  has  been  studied  by 
Hughes  et.al.  [36].  The  theory  includes  frictional  effects  and  allows 
for  geometrically  complicated  contact  surfaces.  Contact  is  not  limited 
to  node  to  node  configurations. 

Hybrid  and  mixed  models  have  also  been  used  to  solve  contact  prob¬ 
lems  in  as  attempt  to  obtain  more  accurate  results.  Kubomura  [37]  de¬ 
veloped  a  method  for  analyzing  contact  problems  which  may  be  either 
frictional  or  frictionless  and  may  involve  extensive  sliding  between  de¬ 
formable  bodies.  It  consists  of  a  hybrid  formulation  including  the  con¬ 
ditions  of  contact,  and  a  scheme  to  determine  the  place  of  contact. 

The  work  of  Hung  and  Saxce  [38]  assumes  that  no  friction  between 
solid  bodies  exists  so  that  the  contact  condition  may  be  expressed  in 
terms  of  displacements  alone.  The  superelement  technique  was  used  to 
reduce  the  number  of  degrees  of  freedom.  Hybrid  triangles  with  12  de¬ 
grees  of  freedom  were  used  to  model  the  deformable  bodies. 

The  mixed  finite  element  method  was  applied  to  two  dimensional 
elastic  contact  problems  by  Tseng  and  Olson  [39] .  Some  of  the  displace¬ 
ment  variables  were  treated  as  natural  boundary  conditions  in  the  con¬ 
tact  region.  The  procedure  includes  sliding  and  adhering  of  contact 
node  pairs. 
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The  representation  of  the  surface  of  contact  has  been  a  subject  of 
special  concern.  A  special  bond  element  was  utilized  by  Schafer  [40]  to 
represent  mathematically  the  disipative  processes  which  occur  during 
contact.  A  relationship  between  the  shear  stress  and  the  frictional  de¬ 
formation  is  assumed  to  be  known,  and  an  element  stiffness  matrix  is  de¬ 
rived  with  the  aid  of  the  principle  of  virtual  work. 

Oden  and  Piries  [41]  propose  nonclassical  friction  laws  in  an  at¬ 
tempt  to  overcome  physical  and  mathematical  difficulties  which  arise 
through  use  for  the  classical  Coulomb  law  of  friction.  Variational 
principles  for  boundary-value  problems  in  elasticity  in  which  such  non¬ 
linear  nonlocal  laws  hold  were  then  developed. 

The  problems  of  elastic  beams  or  plates  supported  on  a  foundation 
have  been  examined.  An  iterative  technique  was  presented  by  Svec  [42] 
for  determination  of  contact  area  between  a  plate  and  its  supporting 
elastic  foundation. 

Mohr  [43]  modelled  an  elastic  supporting  beam  or  plate.  The  method 
is  not  iterative.  A  contact  stiffness  swtrix  is  calculated  to  simulate 
adhesion  to  rigid  surfaces. 

Westbrook  [44]  solved  contact  problems  for  the  elastic  beam.  Vari¬ 
ational  inequalities  and  finite  elements  wore  used  to  obtain  approximate 
solutions  to  bending  problems  for  elastic  beams  in  the  presence  of  a 
rigid  barrier. 

Solutions  of  the  problem  of  wheel  contact  have  been  reported. 
Tielking  and  Schapery  [43]  presented  a  method  of  shell  contact  using 
discrete  Fourier  transform  for  analysis  for  linear  and  certain  types  of 
nonlinear  problems.  The  method  is  used  to  calculated  the  road  contact 


pressure  predicted  by  a  finite  element  toroidal  shell  model  of  a  pneu¬ 
matic  tire. 

A  numerical  method  was  developed  by  Paul  and  Hashemi  [46]  for  the 
determination  of  contact  pressure  which  arises  when  two  elastic  bodies 
with  closely  conforming  non-Hertzian  frictionless  surfaces  are  pressed 
together.  The  method  included  a  technique  for  automatically  altering 
meshes  that  overlay  the  changing  contact  patches.  Problems  of  wheel  and 
rail  contact  were  examined. 

Mathematical  techniques  have  also  been  applied  to  contact  problems. 
Singh  and  Paul  [47]  showed  how  to  solve  anti  formal  non-Hertzian  problems 
using  the  so-called  simply  discretized  method.  The  problems  solved  in¬ 
volved  frictionless  surfaces. 

Conry  and .  Seireg  [48]  used  optimization  techniques  to  develop  a 
mathematical  programming  method  which  utilized  a  simplex-type  algorithm. 
The  technique  is  applied  to  contact  of  beams  of  elastic  foundations. 

The  present  work  is  an  extension  of  the  work  done  by  Kubomura  [37] . 
The  hybrid  stress  model  is  used  to  solve  edge  contact  problems  involving 
thin  multilayer  laminated  composite  plates,  modeled  using  the  newly  de¬ 
veloped  multilayer  elements.  The  model  can  analyze  problems  in  which 
there  may  or  may  not  be  friction  along  the  contact  surface.  There  may 
be  extensive  sliding  between  the  two  contacting  bodies.  "Contact  ele¬ 
ments"  having  unknown  nodal  contact  tractions  are  used  to  model  the  con¬ 
tact  surface. 

In  general,  the  contact  of  an  elastic  body  with  the  edge  of  a  lami¬ 
nated  plate  could  produce  complicated  three-dimensional  effects.  There 
might  include,  for  example,  introduction  of  bending  forces  and  bending 


displacements  even  in  symmetric  laminates,  severe  warping  in  the  zone  of 
contact,  and  a  contact  surface  which  varies  through  the  thickness  of  the 
laminate.  In  principle,  such  a  3-D  contact  analysis  could  be  developed. 
However,  such  an  analysis  was  felt  to  be  impractical  at  this  stage  of 
investigation  because  of  the  excessive  computer  core  and  CPU  require¬ 
ments  needed  to  model  these  effects  (i.e.  a  more  general  multilayer 
plate  element,  and  an  arbitrary  3-D  contact  surface  inspection  and  cor¬ 
rection  scheme).  Therefore,  attention  is  restricted  here  to  thin  symme¬ 
tric  laminates  subject  to  edge  contact  which  does  not  vary  in  the  thick¬ 
ness  direction  and  for  which  the  contactor  moves  only  in  the  plane  of 
the  laminate.  In  effect,  only  inplane  laminate  displacements  will  occur 
even  though  the  bending/stretching  multilayer  element  MQH3T  will  be 
used. 


3.2  HYBRID  STRESS  FUNCTIONAL  FOR  CONTACT 

Variational  Principle  for  the  Contact  Problem 

Let  us  consider  two  elastic  bodies  -  body  A  and  body  B.  Body  A 

and  body  B  are  in  contact  along  the  contact  surface  Sc.  The  boundary 
of  body  A  and  the  boundary  of  body  B  along  the  contact  surface,  Sc  oc¬ 
cupy  the  same  space,  but  the  two  bodies  do  not  overlap. 

As  can  be  seen  in  Figure  9,  each  body  is  under  its  own  system  of 
prescribed  tractions  and  prescribed  displacements.  Tractions  along  the 
contact  surface  are  resolved  into  normal  and  shear  components. 

Using  equilibrium  we  can  delate  these  tractions. 


t/v-  Ts 0  «  <0 


Tr,*+  Tn9  *  <9 


(3.1) 
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Sliding  between  the  two  elastic  bodies  relative  to  one  another  is  al¬ 
lowed.  When  sliding  occurs  the  following  contact  conditions  must  be 
satisfied  for  each  of  the  contacting  bodies  along  the  contact  zone. 

75  -  ± fJ.Tn  (3*2) 

Where  Tn  is  the  normal  traction  on  the  surface  of  contact,  T,  is  the 
tangential  traction  of  the  surface  of  contact,  and  m  is  the  coefficient 
of  friction.  Note  that  (  )*  and  (  )t>  (or  A  and  B)  denote  quantities 
pertaining  to  body  A  and  body  B  respectively. 

The  location  of  the  contact  surface  is  assumed  between  the  two 
elastic  bodies  as  illustrated  in  Figure  10.  At  times  during  contact, 
the  bodies  may  undergo  extensive  sliding  relative  to  one  another.  Other 
studies  assumed  that  nodes  of  one  body  remain  in  contact  with  nodes  of 
the  other.  This  may  not  be  adequate  because  when  two  bodies  come  into 

contact  it  is  likely  that  the  surface  of  one  body  will  move  relative  to 

the  other  body.  In  the  formulation  used  here,  one  contacting  body  is 
allowed  to  slide  relative  to  the  opposing  contacting  body.  This  motion 

of  the  contacting  bodies  causes  unknown  nodal  forces  to  be  applied  to 

the  contacting  nodes.  The  nodal  forces  of  one  contacting  body  are  in 
equilibrium  with  the  nodal  forces  applied  to  the  opposing  body;  the  pro¬ 
cedure  used  here  results  in  coupled  matrix  equations  with  nodal  forces 
within  the  zone  of  contact  as  additional  unknown  quantities. 
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From  the  theory  of  elasticity  certain  conditions  must  be  satisfied 
if  two  elastic  bodies  are  in  contact.  First,  in  the  region  of  contact 
the  tractions  are  equal  in  magnitude  and  opposite  in  direction.  Second, 
the  normal  tractions  are  compressive  and  the  shear  tractions  oppose  the 
direction  of  motion.  Coulomb's  law  of  friction  is  satisfied.  Third, 
there  are  no  gaps  and  no  penetration  in  the  region  of  contact. 

The  conditions  above  are  imposed  upon  stresses  and  displacements  in 
the  region  of  contact  as  follows: 

Displacements 


Stresses 

T»*  +  T»a  =  o  (3-3) 

T*  *  +  Tse  =0 

Ts*  *  ±  TrtA 

±/jlT/ 


where  Xj  are  the  coordinates  before  deformation,  U£  are  the  dis¬ 
placements,  Tn  are  the  normal  tractions  on  the  surface  of  contact, 

Tg  are  the  tangential  tractions  on  the  surface  of  contact,  and  m  is 
the  coefficient  of  friction.  The  convention  of  the  last  of  equations 


•  L  •  «  •  .  *  .  «  V.  »  '  V  L«  •  * 
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(3.3)  is  such  that  tangential  tractions  counteract  relative  movement. 
The  equality  holds  when  there  is  sliding;  otherwise  the  inequality 
holds. 


The  hybrid  stress  functional  [38]  is  based  on  a  modified  complemen¬ 
tary  energy  principle  in  which  interelement  traction  continuity  and  me¬ 
chanical  boundary  conditions  have  been  released  via  the  Lagrange  Multi¬ 
plier  technique,  and  is  given  by 


(3.4) 


where,  is  the  elastic  compliance  tensor,  T£  are  tractions, 

T £  are  prescribed  tractions,  ff£j  are  stresses,  ?n  is  the  volume  of 

the  nth  element.  Sn  is  the  interelement  boundary  of  the  nth  element, 

and  Sa  is  the  portion  of  the  boundary  of  the  nth  element  on  which 
n 

stresses  are  prescribed.  T£“tf£ j v •  where  v ■  are  the  outward  unit 
normal  components,  and  aij,j  “0  from  equilibrium.  Note  the  second 

term  of  equation  (2.1)  is  equivalent  to  the  second  term  of  equation 

(3.4).  This  is  seen  by  applying  the  divergence  theorem  to  the  second 
term  of  equation  (2.1)  and  requiring  that  stresses  satisfy  the  homogene¬ 
ous  equilibrium  equations  (as  required  by  the  hybrid  stress  model). 

The  conditions  of  contact,  equations  (3.3),  are  introduced  into  the 
hybrid  stress  functional  nae,  equation  (3.4),  using  the  Lagrange  mul- 
tiplier  technique,  where  TQ  and  Tg  are  the  Lagrange  multipliers. 

The  hybrid  stress  functional  naec,  which  includes  the  conditions  of 
contact  imposed  upon  the  two  contacting  bodies  within  the  zone  of  con¬ 
tact  is  stated  in  the  form  [38] : 


nryfc 


J7A+t  ^ki  ^  ~J4  g  1  >  ^ 

-s*  *e  C 

JJf  U:ds  ~  /  L(U"  +  K*) ~  (Ur,'+  x"e\ ] 


+  r>  l(UsA+  x**)  -  (Uj+  xs‘)] Jjs 

To  find  the  stationary  value  of  the  functional,  we  take  the  varia¬ 
tion  with  respect  to  displacements,  Uif  stress,  cr^,  and  the  La- 
grange  multipliers,  T#  and  Tn.  The  resulting  Euler  equations  which 


deal  with  the  surface  of  contact  are  given  by: 
Displacements 

-  (un‘  + 

-  (u,s+xsB)  =  o 

Stresses 

(Ty 

(T*  -  T„ft)=d 
(Ji  +■  • 

ft  -  Tsb)  =  0 


(3.6) 


Coulomb's  law  of  friction  is  also  one  of  the  conditions  of  contact 
which  must  be  satisfied.  A  constraint  is  therefore  added  to  satisfy  the 
condition  dealing  with  the  friction  which  occurs  when  the  opposing  bod¬ 
ies  slide  relative  to  one  another;  i.e: 
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(3.7) 


Substituting  equation  (3.7)  into  the  Euler  equations  which  deal  with  the 
surface  of  contact,  equation  (3.6),  yields 


(  Ur>*  f  U  +)  -  ( 0r,a+ 

OJS  +  Xs*)-  (Use  -h  Xi*)-0 


r»  ‘j  «-  o 


(3.8) 


l  ^  j6*)  »  o 


T  4-  y  -r  if 

•*  -  U-Tr. 


7iS^  fJ-Tr)& 

These  conditions  on  the  field  variables  in  the  zone  of  contact  are 
exactly  the  sane  as  the  conditions  of  contact.  Thus,  the  functional 
II^6,  equation  (3.5),  will  satisfy  the  conditions  of  contact  (within 
the  zone  of  contact)  necessary  to  solve  the  problen  of  two  bodies  in 
contact. 


.  *  W  *.«.  '  < 

/' V*  W 


73 


3.2.2  Incremental  Assumed  Hybrid  Stress  Formulation 

Geometric  nonlinearities  in  the  relationship  between  load  and  dis¬ 
placement  complicate  the  general  contact  problem;  i.e.  as  loading  condi¬ 
tions  change,  the  surface  of  contact  also  changes  in  a  nonlinear  fa¬ 
shion.  The  best  method  to  approach  such  a  problem  is  with  an 
incremental  formulation.  As  the  loading  condition  changes  incrementally 
from  state  to  state,  the  place  of  contact  changes. 

It  is  difficult  to  treat  the  contact  surface  as  an  unknown  in  the 
formulation  because  we  are  unable  to  compute  the  necessary  integrals 
over  the  surface  of  contact.  The  method  used  to  solve  the  contact  prob¬ 
lem  involves  assuming  a  contact  surface  and  computing  the  field  vari¬ 
ables  which  render  the  functional  stationary.  It  is  then  determined 
whether  the  field  variables  satisfy  the  conditions  of  contact  within  the 
contact  zone.  If  the  conditions  of  contact  are  not  satisfied,  another 
contact  surface  is  assumed.  This  process  (iteration)  is  repeated  until 
all  conditions  are  satisfied. 

The  nonlinear  an»^ysis  becomes  a  step-by-step  solution  which  has 
been  decomposed  into  a  series  of  linear  incremental  problems  with  ini¬ 
tial  conditions.  For  instance,  a  displacement  u  may  be  expressed  in 
terms  of  an  initial  value,  u®,  and  and  incremental  value,  Au. 

u.°  -hAU  »•*> 


At  the  beginning  of  each  step,  the  initial  value  of  the  field  vari¬ 
ables  are,  in  general,  not  exact.  Hence  they  may  not  satisfy  stress 
equilibrium  conditions  or  the  compatibility  requirements.  Correction 


■  _  PL _ <«  .  < 
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tera*  for  equilibrium  and  compatibility  can  be  included  in  a  finite  ele¬ 
ment  formulation  based  on  an  incremental  variational  principle.  The  in¬ 
cremental  hybrid-stress  functional  for  contact  as  developed  by  Kubomura 
[38]  may  be  stated  as: 


^'3  dV 


AT:  A  U;  <V<T 


-  f  AT;  Au;ds  -  f  77  6  hi  d$  +(  77  AU(dS 


(3.10) 


where : 

is  the  incremental  stress  tensor. 

AT^  are  the  incremental  tractions. 

AU£  are  the  incremental  displacements  on  the  boundary. 

AT£  are  the  incremental  applied  tractions. 

(Tj+ATp  are  the  nodal  contact  forces  within  the  zone  of  contact. 
^  are  the  displacements  on  the  boundary,  dVn. 
is  the  volume  of  the  n^*1  element, 
is  the  boundary  of  the  nth  element. 
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Sa  is  the  portion  of  the  boundary  surface  on  which  tractions  are 
n 

prescribed. 

S„  is  the  portion  of  the  boundary  surface  on  which  displacements 
n 

are  prescribed. 

Sc^  is  the  portion  of  the  boundary  surface  the  two  bodies  have  in 
common,  known  as  the  contact  zone. 

When  using  the  functional  nmc2  the  following  conditions  must  be 
exactly  satisfied: 

u 

Av  =  h  <f,}  ^  (3 

-n  =  ^ 

Au*  AU  Sliri 


where  AU  is  the  incremental  form  of  the  prescribed  displacements . 

The  conditions  of  contact  within  the  zone  of  contact,  previously 
stated  in  equation  (3.3),  are  stated  now  in  an  incremental  form: 


(T»S  +&T*S)~0 
Cn*  -h  An*  )  +  (Tit  ■hATie)=  O 


AT//  i/A. 

It/ A-  AT/) A  ±h  rtf  +  A  7/) 
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(3.12) 


(On**  +  X*  *)  —  (  Urt^-hAOw6  +*  )  =  O 

(U*A  -hAU/  j-Xi,A)  -  (U*B-t-AU<*  +  /a  B)  -  £> 

where 

Tn,  ATn  are  the  initial  and  incremental  tractions,  respectively, 
normal  to  the  contact  surface. 

Ts ,  ATS  are  the  initial  and  incremental  tractions  tangential  to  the 
contact  surface  respectively. 

Un,  AUn  are  the  initial  and  incremental  displacements  normal  to  the 
contact  surface,  respectively. 

Us,  AUg  are  the  initial  and  incremental  displacements  tangential  to 
the  contact  surface,  respectively. 

Xn,  Xg  are  the  normal  and  shear  components  respectively  of  the  coor¬ 
dinates  before  deformation. 

One  important  feature  of  incremental  solutions  is  that  at  the  be¬ 
ginning  of  a  loading  increment  the  solution  is,  in  general,  not  exact, 
hence  it  may  not  satisfy  stress  equilibrium  conditions  and  it  may  also 
violate  the  compatibility  conditions.  The  formulation  of  the  functional 
takes  into  consideration  the  effect  of  initial  equilibrium  im¬ 
balance  and  leads  to  a  system  of  mitrix  equations  for  the  incremental 
solution  with  equilibrium  correction  Violation  of  displacement  compati- 


bility  leads  to  an  incremental  solution  which  also  contains  corrections 
for  the  compatibility  mismatch. 


3.3  MATRIX  EQUATIONS 

The  hybrid  stress  functional  naec,  equation  (3.10),  for  contact 
problems  is  stated  below  in  matrix  form: 


j  i  *  A<Z  oiv 


-  f  ATr£u  d< 

J7-V, 

■+  j  (T  +  AT )T  £0^6 


.  (  TT£dd<; 

Sswt  Ser„ 


Cr-rdf  r  (  Ail*-  Aue)J* 


(T+A?)r  +  .  (U6  +  X*)]U 
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where  (  )i  it  the  transpose  of  (  ) .  The  vectors  end  matrices  are  des¬ 
ignated  by  underscoring  with  a  tilde.  The  functional  nmec,  equation 

(3.13),  includes  equivalent  nodal  forces  for  the  corrections  of  (a)  ini¬ 
tial  material  overlapping,  (b)  separation  and  (c)  initial  force  imba¬ 
lance  at  the  surface  of  contact. 

In  the  hybrid  stress  model  A a  is  interpolated  in  terms  of  a  set  of 
stress  parameters  such  that  the  homogeneous  equilibrium  equations  are 
satisf ied. 


The  boundary  tractions  are  related  to 
Cauchy's  formula): 


(3.14) 

the  same  stress  parameters  (using 


(3.15) 


The  displacements  are  interpolated  in  terms  of  nodal  displacements 
q  and  Aq: 


4"  -  V  if 

$  •  a  4. 


(3.16) 


Coordinates  are  interpolated  in  terms  of  their  nodal  values. 


(3.17) 


UJ±.lU 


■Saibiheaikiaitaiietiitl 


Th*  contact  surface  is  also  discretized  into  finite  elements  called 
contact  elements  having  contact  nodes.  The  contact  tractions  are  inter¬ 
polated  in  the  contact  element  in  terms  of  their  nodal  values,  t  in  the 


form: 


(t  +  &t)  -  ft)  -t 


(3.18) 


We  substitute  equations  (3.14-3.18)  into  equation  (3.13)  to  obtain: 


—  £  I  J  1.6^  Pj  ^  P,  AyS  <d.V 

~j  At  g  tt  A} 

JVr, 

+J  AfiLTrds 

+  [  £ HJ  l r 


T  k/T 


A/  T*  d& 

'V  /V 


(3.19) 


-J  t  u  ti  ( ty4- y*) 

Sr* 

■/  twm  [  y c?  <■**)]  js  I 


Cn 


The  following  matrices  are  then  defined 

-  I  PTSP  <JV 
1-  =  /  £Ti/  otr 

^/-rj  ~ 


=  j  NThJdj 

5dy\ 


1  *  =  f  //r  r 

'"V  V  -"*w 

4v, 


(3.20) 


S  -  / 

ana  i(i/l+&A)-%ipze)]<^ 

Sc, 


™TH  Js 

^<Tj» 

where  Qg  and  Jgn°  *re  stress  equilibrium  checks  on  Sv  and 
Scn»  *a*P«ctively,  and  is  an  initial  mismatch  check  on  Scn. 

Upon  substitution  of  equations  (3.20)  the  hybrid  stress  functional 
for  contact  problems  becomes: 


(3.21) 


V,?  -  sr  U 

+  +■  ^  &« 4  -  4/4r  <%!* 

*-J  ^'irM:s  i^e  -  4y?*r£.-8  4£ 
i- tjf  A&”'  *  4l/r£zs  -  Al‘r  t*e.‘ 


The  element  nodal  displacements  and  contact  element  nodal  tractions 
are  not  independent  from  element  to  element.  The  element  variables,  (), 
can  be  related  to  their  counterpart  assembled  structure  variables,  de¬ 
noted  here  by  ()*,  by  Boolean  matrices,  LnA,  L_B,  and  L_t: 
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Not*  that  th*  element  stress  parameters,  A0,  are  independent  from  ele¬ 
ment  to  element.  It  is  therefore  not  necessary  to  define  a  set  of 
"structure"  stress  paramenters. 

Equations  (3.22)  are  substituted  into  equation  (3.21)  to  yield: 


fi  44 ArtL» A^4~ 

Tz4f:*rtb*£‘  V4' 


-&g rLser&  -£■  &  r/z$/ 

-  i£4T L»*r p-/7 Li r  i-  l* £ , 


(3.23) 


The  stationary  condition  of  ng|Ce  gives: 
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Sine*  A3  are  independent  from  eleaent  to  element  they  can  be  elimi¬ 
nated  on  the  eleaent  level  by  solving  the  first  of  equations  (3.25): 


46*-  tk*''  in1 

6»*  U‘ 


(3.26) 


Substitution  of  equation  (3.26)  into  the  third  and  fourth  of  equations 


(3.25)  yields 
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where 


U  -  4/7  M/"  6/ 


Equation  (3.28)  defines  the  element  stiffness  matrices.  The  summa¬ 
tion  over  n  in  equation  (3.27)  corresponds  to  the  usual  assembly  opera¬ 
tion.  Applying  this  to  equation  (3.27)  and  the  last  of  equations  (3.25) 

gives: 

K**  t  &A* 

<**41**+ f6t+  t*  =  +  a**-  - 

J  (3.29) 

-*£  =  ° 


Equations  (3.29)  may  be  collected  in  matrix  form  in  the  following  man¬ 
ner: 


i<L*A  t-  &+*  -  1 

J>QtS  4.  Q.*6  ) 


(3.30) 


For  the  contacting  bodies,  K£a,  K*B  is  assembled  and  con¬ 
strained,  The  only  portion  of  equation  (3.31)  which  must  be  recomputed 
for  each  iteration  are  the  contact  stiffnesses,  £*A,  _F^B  the  global 
stiffnesses  of  A  and  B  remain  constant.  In  practice  using  the  first  two 
of  equations  (3.29)  to  solve  for  A^A  and  Aij^®  in  terms  of  t^. 

Aq*A  and  Aq*®  are  substituted  into  the  last  of  equations  (3.29). 

Once  t*  was  found,  Ag^A  and  A<i^b  were  found  using  back  substitution. 

The  assembled  contact  stiffnesses  are  very  sparce;  i.e.  only  the  dis¬ 
placement  degrees  of  freedom  which  fall  on  the  contact  surface  will  con¬ 
tribute  to  the  contact  stiffness.  The  sparcity  is  considered  in  order 
to  minimize  the  number  of  multiplications  necessary  to  find  a  solution. 

3. A  THE  CONTACT  ELEMENT  FOR  LAMINATE  EDGE  CONTACT 

As  described  earlier,  the  first  step  to  begin  an  increment  or  iter¬ 
ation  is  to  assume  the  location  of  the  contact  surface,  Sc,  between 
the  two  contacting  bodies  (Figure  10).  The  contact  surface  is  divided 
into  contact  elements  with  contact  nodes  (Figure  11).  It  was  determined 
that  the  elements  should  have  three  nodes  so  that  they  would  be  compati¬ 
ble  with  the  number  of  nodes  along  one  side  of  the  multilayer  laminated 
composite  plate  element  (MQH3T)  which  has  been  described  in  the  previous 
chapter. 

The  normal  and  tangential  tractions  are  interpolated  along  the  con¬ 
tact  surface  in  terms  of  their  values  at  the  three  contact  nodes  using 
the  quadratic  distribution  given  by:  where  |  va. ies  from  -1  to  1,  as 

seen  in  Figure  11,  and  t^  and  tg  are  the  nodal  values  of  the  normal 
and  tangential  components  respectively  of  the  contact  tractions. 


Contact  surface 


3 


Element  (N-l)/2 


Figure  11  Typical  contact  surface  and  description  of 
single  contact  element. 
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-fl-t-ATs  ~  ±  (5* +  (i-X*HsL  +  4  ( 


(3.31) 


The  relationship  described  in  equation  (3.31)  is  stated  in  matrix  form 
by  equation  (3.18),  and  defines  the  matrix  M.  The  contact  element 
stiffness,  F^,  can  then  be  computed  by  integrating  over  the  contact 
surface  using  the  last  of  equations  (3.20). 

In  the  present  formulation,  it  is  not  necessary  for  the  nodes  of 
one  elastic  body  to  be  in  contact  with  the  nodes  of  the  opposing  body 
along  the  contact  surface.  Nor  is  it  necessary  for  the  nodes  of  either 
elastic  body  to  be  in  contact  with  the  contact  nodes.  A  contact  surface 
is  assumed  to  be  fixed  in  the  coordinate  system  but  not  necessarily  to 
the  contacting  bodies.  The  contact  surface  and  the  contact  nodes  are 
assumed  to  be  known  in  order  to  perform  the  integrals  along  the  contact 
surface,  Sc. 

In  application,  the  interval  of  the  integral  is  broken  into  a  sum 
of  integrals  each  of  which  is  over  an  area  in  which  one  contact  element 
is  in  contact  with  one  body  element.  In  Figure  12,  the  interval  of  the 
first  integral  begins  at  node  1  of  contact  element  A  and  corresponds  to 
the  degrees  of  freedom  at  the  3  nodes  of  body  element  a  and  also  to  the 
contact  tractions  at  the  3  nodes  of  contact  element  A.  The  interval  of 
the  first  integral  ends  at  node  3  of  the  body  element  a.  The  second  in¬ 
tegral  begins  at'  the  first  node  of  the  body  element  b  and  ends  at  the 
third  node  of  contact  element  A.  The  third  integral  begins  at  the  first 


node  o£  contact  element  B  and  continues  to  the  third  node  of  body  ele¬ 


ment  b. 

This  process  continues  until  the  end  of  the  contact  surface  is 
reached.  As  the  integrals  are  evaluated  they  are  assembled  to  form  the 
contact  stiffness  F* 
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Fiaura  12  Daaeri ption  of  how  tha  integral  necaasary  to 
coapute  th«  contact  atiffneaa  it  brokan  up  along  tha  contact 
aurfaca. 
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3.5  SOLUTION  TECHNIQUE 

3.5.1  Introduction 

Throughout  the  formulation,  tha  contact  surface  has  been  assumed  to 
be  known.  As  the  load  varies,  a  major  portion  of  the  solution  is  to 
vary  the  contact  surface 

The  overall  strategy  for  solving  the  contact  problem  is  illustrated 
in  Figure  13,  and  is  briefly  summarized  as  follows  (for  the  Nth  incre¬ 
ment  and  the  kth  node) : 

1.  The  Nth  increment  of  external  load  or  prescribed  displacement  is 
applied. 

2.  A  contact  surface  is  assumed  along  with  the  contact  nodes  on  the 
contact  surface.  Each  contact  node  is  assumed  to  be  sliding  or 
not  sliding. 

3.  The  necessary  matrices  are  calculated  and  assembled.  Then  for 
the  ith  iteration  of  the  N^h  load  step,  incremental  displace- 

•  .  /t/ 

ments  and  contact  traction*  1(Tit+ATk)  art  solved 

for  from  the  finite  element  matrix  equation  (3.31), 

4.  Knowing  the  total  displacement  ^'^U^  at  the  end  of  the  previ¬ 
ous  loading  step  (N-l) ,  the  total  displacements  N~^U]t  + 

A*Uk  on  the  boundary  are  cheeked  to  determine  if  the  dis¬ 
placements  satisfy  the  condition  that  there  are  no  gaps  and  no 
penetration  in  the  region  of  contact.  The  contact  tractions 
x(Tk+ATk)  are  checked  to  determine  if  the  normal  tractions 


are  compressive  and  if  Coulomb's  law  of  friction,  equation  (3.2), 
is  satisfied.  If  any  of  these  conditions  are  not  satisfied,  the 
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Figura  13  Suaasary  of  tha  solution  tachniqua  for  solving  tha 
problam  of  contact  batvoan  two  alastic  bodias  using  tha 
hybrid  strass  formulation. 
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location  of  tha  assumed  contact  aurfaea  it  aodifiad  and  stapa  3 
and  4  ara  rapaatad  until  tha  conditiona  of  contact  ara  aatiafiad. 

5.  Once  tha  conditiona  daacribad  abova  hava  baan  aatiafiad  a  conver¬ 
gence  check  ia  aade.  If  tha  aolution  is  not  converged  atapa  3  - 
3  ara  rapaatad  until  tha  aolution  doaa  converge. 

Specifically  tha  conditiona  of  contact  are  liatad  below  with  a 
brief  daacription  of  tha  Banner  in  which  tha  contact  aurfaea  ia  aodifiad 
if  tha  condition  of  contact  ia  not  aatiafiad. 

1.  The  firat  condition  of  contact  atataa  that  tha  nodea  of  one  body 

beyond  the  laat  contacting  node  nay  not  penetrate  tha  other  body. 
If  the  alaatic  bodiaa  in  contact  have  overlapped,  tha  contact 

aurfaea  ia  extended  to  include  tha  contact  nodaa  that  penetrated 
the  opposing  body.  Subsection  3.5.2  will  describe  tha  aathod 
used  to  extend  tha  contact  surface. 

2.  The  second  condition  of  contact  states  that  noraal  tractions 
along  the  surface  of  contact  aust  be  eoapressive.  If  the  noraal 
traction  at  any  two  contact  nodes  within  a  single  contact  eleaent 
are  found  to  be  tensile,  the  contact  eleaent  is  excluded  froa  the 
contact  surface.  Subsection  3.5.3  will  explain  the  aethod  used 
to  release  contact  nodea. 

3.  The  third  condition  of  contact  is  the  relationship  between  noraal 
and  tangential  contact  tractions: 


/?,  *- iTs I £*  JA.  iTy.i'/rf-y,  / 


(3.32) 


If  the  equation  (3.32)  is  not  satisfied  then  the  constraint  in 
equation  (3.33)  is  applied  and  the  contact  node  is  assumed  to  be 
sliding.  Subsection  3. 5. A  will  describe  the  method  to  satisfy 
the  relationship  between  normal  and  tangential  contact  tractions. 


/C  -  ju  !t~  i-  AT*  j 


(3.33) 


Subsection  3.5.5  will  describe  the  method  employed  to  check  for 
convergence. 

3.5.2  Extension  of  the  contact  surface 

As  previously  described  the  first  condition  of  contact  states  that 
nodes  of  neither  body  beyond  the  last  contacting  node  may  penetrate  the 
opposing  body.  Each  node  outside  the  zone  of  contact  along  the  surface 
of  two  elastic  bodies  which  may  come  in  contact  with  the  opposing  body 
is  checked  to  determine  if  it  has  penetrated  the  surface  of  the  opposing 
body. 

The  method  to  determine  if  nodes  of  one  body  have  penetrated  the 
other  is  described  as  follows: 

1.  The  first  step  is  to  check  if  any  nodes  in  body  A  penetrate  body 
B.  The  x-coordinate  of  each  node  on  body  A  is  compared  to  the 


x-coordinates  of  the  nodes  on  the  surface  of  body  B  to  determine 
between  which  two  nodes  on  body  B  the  node  on  body  A  falls.  The 
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y-coordinate  of  the  surface  of  body  B  is  computed  at  the  x-coor- 
dinate  at  which  the  node  of  body  A  is  located.  Thus  it- can  be 

determined  if  the  node  has  penetrated  the  surface. 

2.  The  second  step  is  to  check  if  any  nodes  on  body  B  penetrate  body 
A  in  a  manner  similiar  to  that  which  haa  just  been  explained. 
Once  it  has  been  determined  which  nodes  have  penetrated  the  op¬ 
posing  body  for  both  elastic  bodies,  the  contact  surface  is  ex¬ 
tended  to  include  these  nodes. 

3.  The  third  step  is  to  extend  the  contact  surface.  Figure  14  il¬ 
lustrates  before  and  after  the  solution  for  an  iteration  has  been 
computed.  Node  4  and  node  5  of  each  body  have  penetrated  the  op¬ 
posing  body.  The  contact  surface  is  extended  to  include  nodes  4 
and  5  of  body  B.  The  coordinates  of  the  new  node  4  of  body  B  are 
chosen  by  constructing  a  straight  line  between  the  locations  of 
node  4  of  body  B  for  the  solution  of  the  previous  iteration  and 
the  solution  of  the  present  iteration.  The  new  location  of  node 
4  of  body  B  is  located  along  this  line  midway  between  the  bounda¬ 
ry  of  body  A  and  the  boundary  of  body  B.  The  new  location  of 
node  5  on  body  B  is  found  in  a  similiar  manner. 

The  nodes  of  body  A  which  have  penetrated  body  B  are  then 

placed  on  the  contact  surface  in  the  following  manner.  A 

straight  line  is  drawn  between  the  location  of  node  4  for  the 
previous  solution  and  the  present  solution.  Another  straight 
line  is  drawn  between  the  new  nodes  4  and  5  of  body  B.  The  new 
node  4  of  body  A  is  found  at  the  intersection  of  these  two 

straight  lines  just  mentioned.  The  nodes  of  body  A  which  have 
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penetrated  body  B  beyond  the  last  new  node  on  body  B  are  not  add¬ 
ed,  so  that  the  end  node  on  the  contact  surface  might  always  be  a 
node  on  body  B. 

As  a  result  of  nodes  4  and  5  of  bodies  A  and  B  penetrating 
the  oposing  body,  the  contact  surface  is  extended  to  include 
nodes  4  and  5  of  body  B  and  node  4  of  body  A. 

3.5.3  Release  of  Contact  Modes 

Once  the  solution  has  been  calculated  for  a  given  iteration,  all 
the  contact  nodes  are  checked  to  determine  whether  or  not  the  normal 
traction  is  compressive.  The  contact  elements  used  each  have  three 
nodes  as  described  previously  in  Section  3,4.  The  method  used  to  deter¬ 
mine  whether  or  not  to  exclude  a  contact  element  from  the  contact  sur¬ 
face  was  to  drop  any  element  in  which  the  majority  of  its  nodes  have 
tensile  normal  tractions.  In  other  words,  at  least  two  out  of  three 
nodes  must  have  tensile  normal  contact  tractions  to  exclude  the  contact 
element  from  the  contact  surface. 

3.5.4  Satisfaction  of  the  Relationship  between  Normal  and  Tangential 
Contact  Tractions 

As  the  load  is  applied  incrementally  to  one  or  both  of  the  elastic 
bodies  in  contact,  the  location  of  the  zone  of  contact  will  change.  It 
is  possible  for  the  locations  of  the  contact  nodes  to  also  change  as 
well  as  the  type  of  contact,  i.e.  whether  or  not  the  contacting  bodies 
are  allowed  to  slide  relative  to  one  another  at  any  particular  node. 

To  simplify  programming  efforts  for  this  problem  we  have  chosen  to 


locate  the  contact  nodes  along  the  contact  surface  in  the  same  location 
as  the  nodes  along  the  surface  of  body  B. 


98 


If  no  sliding  occurred  during  the  itn  iteration,  the  location  of 
the  nodes  on  body  A  are  located  by  maintaining  the  ratio  of  AL/L  found 
from  the  converged  solution  for  the  last  load  increment  illustrated  in 
Figure  15  and  stated  as  follows: 


(3.34) 


Therefore  may  be  located.  The  following  test  is  made  at  each 

contact  node. 


!%  U  JU.  /%  1-ATy,  ) 
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If  the  condition  is  not  satisfied  at  a  contact  node  the  following  con¬ 
straint  is  introduced: 


It*  /=  fjL  i%  +  6Ty,] 


(3.36) 


The  state  of  the  contact  node  is  then  changed  to  sliding. 

If  sliding  has  occurred  during  the  ith  iteration,  the  location  of 
the  nodes  of  body  A  are  found  by  projecting  them  perpendicularly  onto 
the  contact  surface.  The  following  test  is  made  at  each  node  on  body  A 
to  determine  whether  of  not  the  node  will  continue  to  slide. 

I A!  At'  I  y  C 


Contact  surface 


a2 


B2  X 


»1  A1 


Converged  solution 
for  the  (N-l)th  load 
increment 


Solution  after  the 

ith  iteration 

of  the  Nth  load  increment 


figure  15  Description  of  ratios  used  to  determine  whether  or  ! 
not  a  node  is  classified  as  sliding  or  not  sliding. 
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where  «  is  some  prescribed  small  number  and  AL',  L'  are  illustrated  in 
Figure  15.  If  this  comparison  is  satisfied,  the  node  will  continue  to 
slide.  If  this  comparison  is  not  satisfied,  the  sliding  condition, 
equation  (3.36),  is  dropped  and  the  nodes  of  body  A  are  relocated  in  the 

manner  previously  described  for  the  state  of  no  sliding. 

Thus  after  each  iteration  the  locations  of  nodes  for  both  of  the 

contacting  bodies  are  found  along  the  zone  of  contact  as  well  as  the 

state,  i.e.  sliding  or  no  sliding  at  each  node. 

3.5.5  Convergence 

Once  the  conditions  of  contact  are  satisfied  a  test  of  convergence 

may  be  made, 

l'»n  ; 

j?  =  /  4 Ujc  -  A  Lk 

/  1*rnZ 

where  AUk  is  the  displacement  at  the  degree  of  freedom.  If  R  is 
less  than  a  prescribed  quantity  the  solution  is  considered  converged. 


-  *  •  '  o 
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3.6  EXAMPLE  PROBLEMS  AND  NUMERICAL  RESULTS 

In  this  section,  the  finite  element  model  and  solution  scheme  de¬ 
veloped  in  this  chapter  is  applied  to  problems  of  contact  between  a  disk 
and  the  edge  of  a  rectangular  plate  (used  to  simulate  a  semi-infinite 
half-space).  This  model  problem  was  selected  for  several  practical  rea¬ 
sons  : 

1.  For  elastic  bodies,  an  analytic  solution  exists  for  the  problem 
of  frictionless  contact  between  a  disk  and  a  semi-infinite  half¬ 
space. 

2.  It  was  determined  that  a  relatively  coarse  mesh  would  be  required 
in  order  to  not  exceed  available  core  storage.  Comparison  with 
an  available  analytic  solution  was  necessary  for  identification 
of  an  adequate  mesh. 

The  thin  multilayer  plate  element,  MQH3T.  developed  in  the  previous 
chapter,  has  been  used  to  model  each  of  the  contacting  bodies.  Atten¬ 
tion  is  restricted  to  symmetric  laminates  (i.e.  no  bending/stretching 
coupling)  subject  to  contact  by  a  cylindrical  isotropic  body.  It  is  as¬ 
sumed  that  contact  varies  only  in  the  x-y  plane  (i.e  independent  of  the 
plate  thickness  direction).  Thus  transverse  displacement  and  cross-sec¬ 
tion  rotations  are  zero.  In  order  to  conserve  computer  core  storage, 
all  transverse  displacements  and  cross-section  rotations  have  been  con¬ 
strained  to  be  zero. 

After  substantial  trial  and  error,  the  mesh  shown  in  Figure  16  was 
selected.  The  semi-infinite  half  space  is  simulated  by  a  finite  plate 
with  the  same  overall  dimensions  of  the  disk  due  to  the  limitation  on 
the  number  of  elements.  There  were  not  enough  elments  available  to  make 
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a  transition  froa  tha  vary  saall  alaaents  naar  tha  contact  zone  to  the 
vary  large  alaaents  far  away  froa  the  contact  zona.  Note  that  core 
storage  of  4096k  was  needed  for  this  aash  even  after  constraining  tha 
transverse  displacaaants  and  cross-section  rotations.  (The  aajor  stor¬ 
age  requireaents  correspond  to  the  asseabled  stiffness  and  the  product 
H~^G  stored  for  each  eleaent.) 

At  this  point  due  to  the  coarse  aesh  which  was  used,  there  is  a 
necessary  lack  of  autoaation  to  find  a  reasonable  solution.  For  each 
increaent  of  load  it  was  necessary  to  Bake  saaple  runs  and  coapare  thea 
to  one  another  in  order  to  find  the  increment  of  applied  displacement 
that  would  bring  the  desired  number  of  contact  elements  into  contact. 
This  was  accoaplished  by  keeping  the  magnitude  of  the  contact  traction 
at  the  last  node  to  coae  into  contact  as  close  to  zero  as  possible.  If 
the  magnitude  of  the  last  contact  traction  was  negative,  the  magnitude 
of  the  increaent  of  applied  displacement  was  manually  increased.  If  the 
magnitude  of  the  last  contact  traction  was  greater  than  the  contact 
traction  at  its  neighboring  node  the  aagnitude  of  the  increment  of  ap¬ 
plied  displaceaent  was  manually  decreased.  For  a  finer  aesh  this  would 
not  be  necessary;  i.e.  load  or  displaceaent  increments  could  be  imposed 
autoaatically. 

The  following  probleas  were  solved: 

1.  Moderate  sliding  frictionless  contact  of  a  deformable  disk  with  a 
deformable  seai-inf inite  half-space. 

2.  Moderate  sliding  frictionless  contact  of  a  nearly  rigid  disk  with 
a  deformable  seai-inf inite  half-space. 


MaiiM 


problem  of 
lf~»paco. 
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3.  Moderate  sliding  frietionless  contact  of  a  nearly  rigid  disk  with 
a  symmetric  multilayer  composite  plate. 

In  the  following  subsections,  the  solution  of  each  problem  is  dis¬ 
cussed  in  detail  and  compared  with  independent  solutions  when  available. 

3.6.1  Contact  between  an  Elastic  Disk  and  an  Elastic  Half-space 

The  problem  of  frictionless  contact  between  an  elastic  disk,  and 
an  elastic  half-space  has  been  solved  by  Hertz  [49] .  The  problem  is 
also  solved  in  this  section  by  the  finite  element  method  and  results  are 
compared  with  the  Hertz  solution. 

The  mesh  and  boundary  conditions  are  shown  in  Figure  16.  The  ma¬ 
terial  was  isotropic  and  the  material  properties  for  both  the  disk  and 
the  half-space  are  E“lxl0^  psi.  and  t».25.  Quantities  of  interest  are 
the  contact  tractions,  (T+AT)  and  the  distributions  of  cy  in  the 
half-space  near  the  boundary  (y*0)  and  near  the  plane  of  symmetry  (x-0) . 

The  solution  is  presented  for  the  case  of  prescribed  displacements 
applied  to  the  top  of  the  disk.  The  prescribed  displacements  were  ap¬ 
plied  in  three  increments.  For  the  first  increment  equal  to  -.9698  in. 
two  elements  are  in  contact,  the  second  increment  equal  to  -.148  in.  has 
four  elements  in  contact  and  the  third  increment  equal  to  -.0942  in.  has 
five  elements  in  contact. 

Figure  17  illustrates  the  distribution  of  the  contact  traction, 
(T+AT)  along  the  zone  of  contact  as  they  compare  with  the  Hertz  solu¬ 
tion.  The  contact  tractions  calculated  only  at  2  point  Gauss  stations 
of  the  contact  element  are  shown  in  Figure  18.  The  2x2  Gauss  stations 
were  shown  in  the  previous  chapter  to  be  "optimal  sampling  points”.  The 
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magnitudes  of  the  contact  tractions  calculated  at  these  points  agree 
very  closely  with  the  Hertz  solution. 

The  contact  surface  traction  are  calculated  by  a  second  method  us¬ 
ing  the  stress  distributions  and  the  stress  coefficients  A@'s  assumed 
within  each  plate  element.  The  distribution  of  the  contact  tractions 
are  compared  with  the  Hertz  solution  in  Figure  19.  Once  again  the  dis¬ 
tribution  of  contact  traction  agrees  well  with  the  Hertz  solution. 

Stresses  in  the  half-space  near  the  surface  at  y”0  are  shown  in 
Figure  20.  They  are  calculated  from  the  plate  elements  at  the  2x2  Gauss 
stations  nearest  the  y”0  boundary.  The  stresses  within  the  half-space 
near  the  surface  are  only  slightly  smaller  than  the  tractions  on  the 
surface  as  would  be  expected.  The  stress  tends  toward  zero  as  x  in¬ 
creases.  Stresses  in  the  half-space  near  the  plane  of  symmetry  are 
shown  in  Figure  21.  They  are  also  calculated  at  the  2x2  Gauss  stations 
in  the  plate  elements.  The  stress  is  maximum  at  the  surface  and  de¬ 
creases  to  a  minimum  value  as  y  decreases.  The  accurate  prediction  of 
stress  within  the  half-space  for  this  problem  shows  the  mesh  is  fine 
enough  so  that  one  may  expect  good  results  for  problems  which  have  no 
analytic  solution. 
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Figure  17  Contact  traction  (T+AT)  distribution  vs  x  for 
contact  batvaan  a  disk  and  a  half-space  having  identical 
aaterial  properties.  The  displaeeaents  were  applied  in  3 
increaents. 
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Figure  19  Distribution  of  ths  inplan*  normal  stress,  Oy  in 
tha  half-spac*  calculated  along  the  surface  y“0  vs  x.  The 
displacements  were  applied  in  three  increments. 
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Figure  21  Distribution  of  the  inplane  normal  stress,  oy  in 
the  half-space  calculated  through  the  2x2  Gauss  stations 
nearest  the  plane  of  symmetry  vs  y.  Three  increments  of 
applied  displacement  were  used.  Edisk^haif-space*1  • 
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3.6.2  Contact  between  a  Nearly  Rigid  Disk  and  an  Elastic  Half-apace 

The  problea  of  frictionless  contact  between  a  nearly  rigid  disk, 
under  an  applied  load,'  and  an  elastic  half-space  can  be  solved  using  the 
saae  Hertz  solution  described  in  the  last  subsection. 

The  aesh  and  boundary  conditions  are  shown  in  Figure  16.  The  aa~ 
terials  of  the  disk  and  elastic  half-space  are  isotropic  and  the  aateri- 
al  properties  are  Ed£lk“l.xlO^  psi.,  Ehalf-plane“1,x1^  P***  Uld 

The  solution  presented  is  for  the  case  of 
prescribed  displaceaents  at  the  top  of  the  disk  in  three  increments. 
The  aagnitudes  of  the  applied  displaceaents  were  the  same  as  those  used 
in  the  last  exaaple  problea.  They  are  -.0698  in.,  -.148  in.,  and  -.0942 
in.  for  a  total  of  -.312  in.  Quantities  of  interest  are  the  contact 
tractions,  (T+AT)  and  the  distributions  of  in  the  half-plane  near 
the  boundary  (y*0)  and  near  the  center  (x“0) . 

Figure  22  shows  the  distribution  of  the  contact  tractions  calculat¬ 
ed  and  the  distribution  given  by  the  Hertz  solution.  The  pointwise  dis¬ 
tribution  of  the  contact  tractions  does  not  agree  as  well  as  the  previ¬ 
ous  example  with  the  Hertz  solution  although  the  solution  calculated  at 
the  2  point  Gauss  stations  shown  in  Figure  23  is  just  as  accurate  as  in 
the  previous  example.  This  serves  to  point  out  the  importance  of  using 
these  "optimal  sampling  points”  to  calculate  contact  tractions. 

The  contact  surface  tractions  are  calculated  by  a  second  method  us¬ 
ing  the  stress  distribution  and  stress  coefficients  A$'s  assumed  within 
each  plate  element.  The  distribution  of  the  contact  tractions  is  com¬ 
pared  with  the  Hertz  solution  in  Figure  24.  This  time  the  distribution 
of  the  contact  tractions  agrees  well  with  the  Hertz  solution. 
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Figure  23  Contact  traction*,  ,  (T+&T)  distribution  v»  x 
calculated  at  the  2  pt.  Gauss  stations  along  the  contact 
surface.  The  displacements  were  applied  in  3  increments. 
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Figure  24  Distribution  of  the  inplane  normal  stress,  Oy  in 
the  half-space  calculated  along  the  surface  y“0  vs  x.  The 


displacements  were  applied  in  three  increments. 
Edisk/Ehalf-space“1000* 
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Stresses  in  the  half-plane  near  the  surface  at  y“0  are  shown  in 
Figure  25.  They  are  calculated  from  the  stress  coefficients,  A3  at  the 
2x2  Gauss  stations  nearest  the  y~0  boundary.  The  stresses  within  the 
half-space  near  the  surface  are  only  slightly  smaller  that  the  tractions 
on  the  surface  as  would  be  expected.  The  stress  tends  toward  zero  as  x 
increases.  Stress  in  the  half-plane  near  the  plane  of  symmetry  are 
shown  in  Figure  26.  They  are  also  calculated  from  stress  coefficients 
at  the  2x2  Gauss  stations.  The  stress  is  maximum  near  the  surface  and 
decreases  to  a  minimum  value  as  y  decreases. 


calculated  through 
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3.6.3  Contact  between  a  Nearly  Rigid  Disk  and  a  Symmetric  Multilayer 
Composite  Plate 

The  problem  of  frictionless  contact  between  a  nearly  rigid  disk, 
under  an  applied  displacement,  and  a  symmetric  thin  multilayer  composite 
plate  is  solved  in  this  subsection. 

The  mesh  and  boundary  conditions  are  shown  in  Figure  16.  The  ma- 
terial  of  the  disk  is  isotropic  and  the  material  properties  are 
E*2.xl0^®  psi.  and  v«.21.  The  plate  is  a  4-layer  90/0/0/90  cross-ply 
(  6  measured  from  the  x-axis).  The  thickness,  h  of  the  plate  and  the 
disk  is  equal  to  .02  in.  Each  layer  is  of  equal  thickness  h/4,  and  lay¬ 
er  material  properties  are  En”20xl0®  psi.,  E22**2.1xl0^  psi., 
*'12”*'23“-^*  *06  85x10^.  Magnitudes  of  the  three  incre¬ 

ments  of  applied  displacement  are  -.05313  in.,  -.12009  in.,  and  -.07802 
in.  (equal  to  a  total  applied  displacement  of  -.25124  in.) 

Quantities  of  interest  are  the  contact  tractions,  (T+AT)  and  the 
distributions  of  ay  in  the  plate  near  the  boundary  (y~0)  and  near  the 
plane  of  symmetry  (x«0) . 

Figure  27  illustrates  the  distribution  of  the  contact  tractions 
calculated  using  the  finite  element  program.  The  solution  calculated  at 
the  2  point  Gauss  stations  is  shown  in  Figure  28.  Using  the  experience 
gained  in  the  last  two  examples,  these  are  the  best  locations  at  which 
to  evaluate  the  distribution  of  contact  tractions. 

A  weighted  sum  of  the  stress,  Oy  in  all  the  layers  of  the  multi¬ 
layer  composite  plate  have  been  calculated  along  y"*0  using  the  stress 
distibutions  assumed  within  each  element.  The  distribution  of  stress 
shown  in  Figure  29  is  very  close  to  the  distribution  of  contact  traction 
calculated  at  the  2  point  Gauss  stations  in  Figure  28,  although  the 
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Figure  27  Contact  surface  tractions,  (T+AT)  distribution  vs 
x  for  contact  between  a  nearly  rigid  disk  and  a  symmetric 
multilayer  composite  plate.  Three  increments  of  applied 
displacement  were  used. 
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Figure  28  Contact  traction*,  (T+AT)  distribution  v*  x 
calculated  at  tha  2  pt.  Gauss  stations  along  tha  contact 
surface  for  tha  problem  of  contact  between  a  nearly  rigid 
disk  and  a  symmetric  multilayer  composite  plate.  The 
displaements  were  applied  in  3  increments. 


strata  distribution  within  tha  alamant  is  discontinuous  at  interelement 
boundaries .  Tha  stress,  Cy  calculated  along  y“0  using  the  element 
stress  distribution  in  each  of  layers  1  and  2  are  shown  in  Figures  30 
and  31  respectively. 

A  weighted  sum  of  the  stress  in  all  the  layers  of  the  multilayer 
composite  plate  near  the  surface  at  y”0  are  shown  in  Figure  32.  They 
are  calculated  from  the  stress  coefficients,  A0  at  the  2x2  Gauss  sta¬ 
tions  nearest  the  y”0  boundary.  The  stress,  Oy  calculated  at  the  2x2 

Gauss  stations  nearest  the  y“0  boundary  are  shown  for  layers  1  and  2 
separately  in  Figures  33  and  34.  The  stress  tends  toward  zero  as  x  in¬ 
creases.  Stresses  in  the  multilayer  composite  plate  near  the  plane  of 
symmetry  are  shown  in  Figure  33.  They  are  also  calculated  from  stress 
coefficients  at  the  2x2  Gauss  stations.  The  stress,  <7y  calculated  at 
the  2x2  Gauss  stations  nearest  the  plane  of  symmetry  are  shorn  for  lay¬ 
ers  1  and  2  separately  in  Figures  36  and  37.  The  stress  is  maximum  near 
the  surface  and  decreases  to  a  minimum  value  as  y  decreases. 
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Figure  30  Contact  tractions,  (T+AT)  distribution  in  layar  1 
vs  x  calculated  at  the  2  pt.  Gauss  stations  along  the 
contact  surface  for  the  problem  of  contact  between  a  nearly 
rigid  disk  and  a  symmetric  multilayer  composite  plate.  The 
displacements  were  applied  in  3  increments. 


#2  Pt.  Gauss  stations 


x  (in.) 


Figure  31  Contact  traction,  (T+AT)  distribution  in  layer  2 
vs  x  calculated  at  the  2  pt.  Gauss  stations  along  the 
contact  surface  for  the  problem  of  contact  between  a  nearly 
rigid  disk  and  a  symmetric  multilayer  plate.  The 
displacements  were  applied  in  3  increments. 


Figure  32  Distribution  of  the  inplane  normal  stress,  in 
the  multilayer  composite  plate  calculated  through  the  2x2 
Gauss  stations  nearest  the  contact  surface  vs  x.  Three 
increments  of  applied  displacement  were  used. 
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£2  Distribution  or  the  inplane  normal  stress,  < 
2  of  the  multilayer  composite  plate  calculated 
2  Gauss  stations  nearest  the  contact  surface 
increments  of  displacement  were  applied. 
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Figure  3^  Distribution  of  the  inplane  normal  stress,  oy  in 
the  multilayer  composite  plate  calculated  through  the  2x 
Gauss  stations  nearest  to  the  plane  of  symmetry  vs  y.  Thre 
increments  of  applied  displacement  were  used. 
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Chapter  IV 

CONCLUSIONS  AND  SUGGESTIONS  FOR  FURTHER  RESEARCH 

4.1  CONCLUSIONS 

The  purpose  o£  this  study  has  been  to  solve  the  problem  of  edge 
contact  between  a  nearly  rigid  disk  and  a  multilayer  laminated  composite 
plate.  A  suitable  plate  element  has  been  formultated  for  this  purpose 
as  well  as  the  method  to  solvs  the  contact  portion  of  the  problem. 

The  hybrid  stress  formulation  for  an  isoparametric  thin  to  moder¬ 
ately  thick  multilayer  laminated  composite  plate  element  has  been  pre¬ 
sented.  The  displacement  behavior  was  characterized  by  laminate  refer-- 
ence  plane  inplane  and  transverse  displacements  and  laminate  non-normal 
cross-section  rotations;  as  a  result,  the  number  of  degrees  of  freedom 
was  independent  of  the  number  of  layers.  All  components  of  stress  were 
included. 

A  thin  to  moderately  thick  multilayer  laminated  composite  plate  el¬ 
ement  (MQH3)  was  developed.  Due  to  the  excessive  CPU  time  necessary  to 
calculate  the  stiffness  of  a  single  element,  several  unsuccessful 
schemes  were  developed  to  decrease  the  CPU  time  .  An  element  formula¬ 
tion  for  thin  laminates  was  developed  as  a  reduction  of  the  moderately- 
thick  formulation  by  neglecting  the  contribution  of  transverse  shear 
stresses  and  normal  stress  to  the  internal  complementary  energy.  Using 
the  reduced  formulation,  an  algorithm  was  developed  which  decreased  the 
CPU  time  necessary  to  calculate  the  stiffness  by  a  significant  amount. 
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CPU  tine  for  the  reduced  element  is  essentially  independent  of  the  num¬ 
ber  of  layers. 

For  thin  laminates  the  performance  of  the  two  elements  is  essen¬ 
tially  identical.  For  thickness  ratios  h/L>.01,  the  more  general  ele¬ 
ment  is  necessary  to  produce  accurate  results,  but  the  reduced  element 
results  in  a  large  reduction  in  CPU  time.  The  MQH3T  element  retains  a 
high  degree  of  accuracy  for  thin  plates  having  a  thickness  ratio  less 
than  h/L'-O.Ol.  For  the  case  of  an  angle-ply  plate  which  is  simply  sup¬ 
ported  on  all  edges  the  results  are  very  good  and  6  has  little  effect  on 
the  accuracy  of  the  predictions.  The  poor  results  displayed  for  small  6 
in  the  case  of  a  angle-ply  plate  which  has  been  clamped  on  all  edges  are 
thought  to  be  due  to  a  type  of  "material"  aspect  ratio  effect. 

The  hybrid  stress  formulation  for  the  analysis  of  laminate  edge 
contact  problems  as  well  as  the  procedure  to  locate  the  surface  of  con¬ 
tact  has  been  presented.  A  contact  surface  was  assumed  between  two  con¬ 
tacting  bodies  and  was  divided  into  contact  elements.  Attention  was  re¬ 
stricted  to  cases  of  inplane  displacements  and  contact  which  was  uniform 
in  the  through  thickness  direction. 

Examples  considered  were  the  problem  of  contact  between  a  disk  and 
a  half-plane  involving  both  elastic  and  nearly  rigid  disks,  and  isotrop¬ 
ic  and  multilayer  composite  plates.  Results  of  moderate  sliding  contact 
are  in  very  good  agreement  with  Hertz  solutions  and  can  be  considered 
good  enough  to  obtain  stresses  accurately.  The  problem  of  moderate 
sliding  frictionless  contact  of  an  elastic  disk  with  an  elastic  semi-in¬ 
finite  half-space  agreed  very  well  with  the  Hertz  distribution  of  con¬ 
tact  stress.  The  best  distributions  of  contact  tractions  were  obtained 
at  the  2  point  Gauss  stations.  The  results  obtained  from  the  example 


problems  involving  ths  semi-infinite  half-space  suggest  that  the  mesh 
was  fine  enough  to  obtain  reasonable  results  for  the  problem  of  edge 
contact  between  a  disk  and  a  thin  multilayer  composite  plate. 


A. 2  SUGGESTIONS  FOR  FURTHER  RESEARCH 

Suggestion  for  further  research  include: 

1.  Efforts  should  be  made  to  decrease  the  storage  requirements  i.e. 
out-of-core  solvers  and  use  of  disk  storage,  so  that  finer  meshes 
could  be  used  to  solve  more  complicated  problems  such  a  curved 
boundaries . 

2.  Extension  of  the  present  program  to  include  transverse  displace¬ 
ments  and  rotations. 

3.  Use  of  more  sophisticated  multilayer  elements  near  the  contact 
surface  which  would  include  the  possibility  of  cross-sectional 

i 

warping  at  the  edge  of  the  laminate. 

Inclusion  of  material  and  geometric  nonlinearities. 
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